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ABSTRACT 

In this paper, a third in a series we continue the construction and the analysis 

of essentially non-oscillatory shock capturing methods for the approximation of 

hyperbolic conservation laws. We present an hierarchy of uniformly high order 

accurate schemes which generalizes Godunov’s scheme and its second order accurate 

MUSCL exten sion to arbitrary order of accuracy. The design involves an essentially 

non-oscillatory piecewise polynomial reconstruction of the solution from its cell 

averages, time evolution through an approximate solution of the resulting initial 
value problem, and averaging of this approximate solution over each cell. The 

reconstruction algorithm is derived from a new interpolation technique that when 
applied to piecewise smooth data gives high-order accuracy whenever the function 
is smooth but avoids a Gibbs phenomenon at discontinuities. Unlike standard 
finite difference methods this procedure uses an adaptive stencil of grid points and 
consequently the resulting schemes are highly nonlinear. 
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1. introduction. In this paper, the third in a series we continue to study 
the use of essentially nonoscillatory, uniformly high-order accurate schemes for the 
numerical approximation of weak solutions of hyperbolic systems of conservation 


(1.1a) 


ut 4- /(u) x = 0 


(1.1b) 


u(z,0) = u 0 (i). 


Here u = (u t> ... , u m ) T is a state vector and /(u), the flux, is a vector valued 
function of m components. The system is hyperbolic in the sense that the m x m 


Jacobian matrix 


has m real eigenvalues 


A(u) = df/du 


a t (u) < a 2 (u) < < a m (u) 


and a complete set of m linearly independent right-eigenvectors {rjfc(u)}£* =1 . We 
denote by {/fc(«)}/k*=i the left-eigenvectors of A(u ) and assume that l,r k = t k- 

We assume that the initial value problem (IVP) (1.1) (embedded in an appro- 
priate setting which includes entropy considerations) is well-posed in the sense of 
Cauchy and that its weak solutions are generically piecewise smooth. We denote 
its evolution operator by E(t), i.e. 


u( ,0 = E{t) Uq. 


Let tS(r) denote the sliding average of u/(z) 


(1.3a) 


i 

)(z) = - / u>(z + y )dy = {A h w){x) 
* J -hi? 
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We note that d> is smoother than «, by one derivative, and that at points of 


smoothness 


(1.3b) 


U>(z) = w(x) + 0(h 2 ). 


The sliding average in z of a weak solution of (1.1), u(x, t), satisfies 

dt u ( x,t ) + £^( u ( x + V 2 »0) ~ /(«(* - A/2,i))] = 0. 
Integrating this relation from t to t + r, we get 

(1.5a) «(*, t + r) = «(*, t) - A[/(x + A/2, t; «) - /(, - * /2 , u) | 


where A = r/h and 


(1.5b) 


. i r 

f(x,t;w)=- f(w(x, t + ij))drj. 
7 Jo 


Let {I, x [*«,f« + ij/ where I 3 = [x } _^ Xj+ ^ Xa = ah> t? 
partition of R x R + . Writing relation (1.5) at x = x J} t = t n we get 


= nr, be a 


(1.6a) 


u j + U 1 A (/( x j+A ( < n ;u) - /(z y _i,f n ;u)). 


(1.6b) 


*/ = U(ly,t n ) = U(l,t n )d2 


is the “cell-average” of u at time t n . 

In this paper we describe a class of numerical schemes that generalises 
Godunov’s scheme [S| and its second order extensions ((22), |4|, |15|) to any finite 
order of accuracy. These schemes can be written in standard conservation form 

d.7a) 


Here £ft(r) denotes the numerical solution operator and the numerical flux, 

denotes a function of 2k variables 


(1.7b) 


= f( v ?-k+i>-,v"+ k h 


which is consistent with the flux f (u) in (11) in the sense that /(u.u, = f( u ). 
We design these schemes so that the conservation form (1.7a) will approximate (1.5) 
to high order of accuracy. Setting v* = ti? in (1.7) and comparing it to (1.6) we 
see that if the numerical flux (1.7b) can be expanded as 


( L8a ) />-$ = f(x 3 ~ ^ * 0(h r * 1 ) 

then 

- X[d(x } ^x) - d{z,_i) h r ~ 0{h r ~ l ). 

This shows that if the numerical flux Jj + x satisfies (1.8a) then the truncation error 
in the sense of cell averages is 


(1.8b) 


ti(i,.f n * r) - 'E h [r)a( .t n ) j = \\d(z k ) - d(x : _x )'h' * Oih" 


which is 0{h T * 1 ) where d(z) is Lipschitz continuous. 

When /( u) is a nonlinear function of u. the approximation of /'r,_ * . t„. u. t: 
O(h') requires knowledge of pointwise values of the solution to the same order of 
accuracy. In order to design a numerical flux that satisfies (1.8a). we must extract 
high order accurate pcintwise information from :he given {v*}, which are approxi- 
mations to {u* }. the cell averages (1.6b) of the solution. Solving this reconstruction 
problem to arbitrarily high-order of accuracy r, without introducing O 1 Gibbs- 
like spurious oscillations at points of discontinuity, is the most important step in 
the design of our new schemes. 


G!ve n a,- - ®(xy), cell averages of a piecewise smooth function «,(*), W e con 
struct R( r;®), a piecewise polynomial function of x of unifomt polynomial degre, 


( r — 1) that satisfies: 

(i) At all points x for which there is a neighborhood where w is smooth 

R(z;&) = tu(x) + e(x)h r 4 - 0{h r+l ). 

(ii) conservation in the sense of 


(1.9b) 


R(xj\ w) = tDj 


here R denotes the sliding average (1.3) of R. 

(iii) ^ is essentially non-oscillatory 


(1.9c) 


TV (R( ; u>)) < TV(w) + 0(h r ), 


where TV denotes total-variation in x. 

The inequality (1.9c) implies that the reconstruction R is essentially non- 

oscillatory in the sense that it does not have a Gibbs-like phenomenon of generating 

0(1) spurious oscillations at points of discontinuity that are proportional to the size 

of the jump there. In ( 16|, |„|, (,7) we describe *(x;in) in the scalar case. We show 

there that R may occasionally produce 0(K) spurious oscillations which are on the 

level of the truncation error. These small spurious oscillations may occur only in 

the smooth part of „ and they usually disappear once «(,) is adequately resolved 

on the computational mesh. Forsake of completeness we review this reconstruction 

algorithm in section 3; we shall extend it to vector functions ui(x) in Section 5 of 
this paper. 

Using the reconstruction (1.9) we can express the abstract form of our new 
schemes by 


£/»(r) u> = A h E(t) R{ -w). 


( 1 . 10 ) 
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Here A h is the cell-averaging operator on the RHS of (1.3); E(t) is the exact evolu- 
tion operator (1.2) and u> is any piecewise smooth function of x. These schemes are 
a generalization of Godunov’s scheme and its second-order extensions in the sense 
that (1.10) with the first order accurate piecewise constant reconstruction 

O-H) R(x; tD) = tDy for x _i < x < x, . i 

J 7 2 

is exactly Godunov’s schemes ([5]); (1.10) with the second order accurate piecewise 
linear reconstruction 

(1.12a) R[x\ w) = Wj 1 Sj(x - xy) for x < x < x, . i 

J 2 J ’ 2 

such that 

( L12b ) s j — w x(xj) + O(h), 

is the abstract form of the second-order extensions to Godunov’s scheme described 
in [22], [4] and (15). 

We recall that the evolution operator E(t ) is monotone in the scalar case. Since 
Ah, the cell-averaging operator is also monotone we see that in the scalar case 

(1.13a) TV(E h (r)w) = TV (A h E(r) R( -w )) < TV(R( ;w)). 

If tZ> in (1.13a) is the sliding average of a piecewise smooth function w(x), it 
follows then from (1.9c) that 

(l-13b) TV(E h (r)w) < TV(w) + 0(h r ). 

This shows that the schemes (1.10) in the scalar case are essentially non- 
oscillatory in exactly tiae same way as the reconstruction: They do not have a 
Gibbs-like phenomenon at discontinuities, yet they may occasionally produce small 
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spurious oscillations on the level 0{h r ) of the trunction error (see Remark 1.3 at 
the end of this section). 

(1.10) is the abstract operator expression of a scheme in the conservation form 
(1.7). Although the scheme generates discrete values «*, which are r-th order 
accurate approximations to the cell-averages u?, its operation involves a globally 
defined pointwise approximation to u(z ,f) of the same order of accuracy, which we 
denote by v h (x,t). The latter is defined for all z in the time-strips t n < t < t n+ u 
with a possible discontinuity at {tk}\ we shall use the standard notation Vh{x,t K ±0) 
to distinguish between the two possibly different values. 

We define v h (z,f) via the following algorithmic description of the scheme (1.10). 

We start by setting 

v° = Uo(Zj) 

where u 0 is the given initial datum (1.1b), and u 0 is its sliding average (1.3a). 
Having defined v n = (v”), approximation to {u^} in (1.6b), we proceed tc evaluate 
v n+1 by the following three steps: 

(i) Reconstruction : Define 


(1.14a) v h {x,t n +0) = R{x\v n ), 

Note that v h {x,t n +0) is a pointwise approximation to u{x,t n ). 

(ii) Solution in the small : For t n < t < t n + r = t n + i i define 

(1.14b) Vh{ ;<) = E{t - t n ) v h { -,t n + 0). 


(iii) Cell-averaging: Close the time-loop oi the algorithm by defining 


(1.14c) 


v n } JrX = M *>;<»+ 1 - o ) = 



X,t n + 1 - 0 )dx. 


I 
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We note that v*, being an exact solution of (1.1) in t n < t < *n+ 1 > satisfies 
(1.5) in this strip. Using the conservation property (1.9b) of the reconstruction in 
(1.14a), i.e. 

(1.15) v h (zj,t n + 0) = v?, 

we get from (1.5) that the scheme (1.10), (1.14) satisfies the conservation form 
(1.16a) v ; +1 =^-A(/ J+t -/ i _x) 

with the numerical flux 

(1.16b) f J+ i = f(x J+ j,tn;v h ) = ^ j f{v h {x : + i,tn + V))dri. 

We turn now to examine the local truncation error of the scheme. For this 
purpose we consider a si ogle application of (1.14) starting with v " = , the exact 

cell-averages of the solution. It follows from (1.9a) and (1.14a) that 

(1.17a) Vh{x,t n +0) = v(x,t„) + e(z)h r + 0(/t r+1 ). 

The definition (1.14b) and our assumption of the well-posedness of the IVP (1.1) 
imply that 

(1.17b) v h (z,t) = u(zj) + 0(h r ) for t„ < t < 

This in turn * 's that the numerical flux (1.16b) of the scheme satisfies (1.8a). 
i.e. 

(1.17c) i + d{Zj+i)h r +0{h r * 1 ). 

Clearly non-smoothness of d(z) in (1.17c) can result only from non-smoothness 
of the coefficient «(x) in (1.17a). It follows then from (1.8b) that away from points 
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of discontinuity and points at which e(x) fails to be Lipschitz continuous, the local 

truncation error in the sense of cell-averages is 0(h +1 ). 

Let u(x, t) be a smooth solution of (1.1) and let us suppose that as h -» 0, r = 
O(h), the numerical approximation converges pointwise to u(x,t). If e(x) is globally 
Lipschitz cntinuous then the local truncation error in the sense of cell averages is 
globally 0(V +1 ). At time t, after performing N = t/r time-steps, we expect the 

cumulative error to be 0(/i r ), i.e. 


(1.18a) 


Vj* = u(x_,,fjv) + 0(h r ) 


In this case we see from (1.9a) that 


(1.18b) 


v h {x,t N + 0) = R{x\v n ) = u{x,t N ) + 0{h r ). 


Thus at the end of the computation we have two sets of output data at our disposal, 
(i) Discrete values {v^} that approximate {u(xj,t N )} to 0[h r ). (n) A piecewise 
polynomial function of x, R{x\v N ), that approximates u{x,t N ) to 0{h r ). 


REMARK: (1.1). Note that (1.8) is quite different from the truncation error in 
a pointwise sense which is used in formulating Lax- Wendroff- type schemes [20], 
[21]. There we take v? = u(xy,t n ) and require v" +1 = u(xy,t n +i) + 0(/i r+l ). To 
accomplish that we need a numerical flux that satisfies 




fc= 1 


We shall see in the following that condition (1.8a) for the accuracy in a cell-average 
sense is more manageable in many respects 


•1 


REMARK (1.2): When e(x) fails to be Lipschitz continuous at a point, the local 
truncation error (1.8b) is only 0{h r ). In the MUSCL-type schemes [22], [4] this 
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happens at local extremum points' in higher-order accurate schemes this may occur 
at roots of higher derivatives of u (see [15], [11]). Due to local accumulation we 
expect the pointwise error at time t , after N = tjr time-steps, to be only 0(h r ~ 1 ) 
at such points. Away from these points we expect the pointwise cumulative error to 
remain 0(h r ). Consequently the scheme is (r- l)-th order accurate in the maximum 
norm. Because of the non-oscillatory nature of the schemes, we expect the number 
of points where e(x) fails to be Lipschitz-continuous to remain bounded as h —* 0. 
In this case the Li-norm of the cumulative error is 0(h r ). To distinguish between 
schemes that are r-th order accurate in the usual pointwise sense, and those that 
are r-th order accurate in the ii-norm but only (r- l)-th accurate in the maximum 
norm, we shall use “r-th order accurate” for the latter, thus qualifying the difference 
by the use of quotation marks. 

REMARK ( 1.3) : It is well known that if the total variation of the numerical ap- 
proximation is uniformly bounded, i.e. 

(1.19) TV(v h (- y t))<C TV (u 0 ) 

where the constant C is independent of h for 0 < t < T, then any refinement 
sequence h — ► 0, r = 0(h) has a subsequence that converges in L\ cc to a weak 
solution of (1.1). Therefore uniform boundedness of the total variation is an appro- 
priate sense of stability for numerical approximations to discontinuous solutions of 
(1.1); see [9], [10] and the references cited there. 

Inequality (1.13) shows that the total variation of our new schemes is dominated 
by that of reconstruction step 


( 1 . 20 ) 


JV(</" +1 ) < TV{R( 
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When R is the piecewise-constant function (1.11) or the piecewise-linear function 
(1.12) (where the slope 3j is that of the MUSCL scheme) then 

(1.21a) TV{Ri-,v))<TV{v) 

for any function v of bounded total variation. Consequently Godunov’s scheme and 
the M T JSCL scheme are total variation diminishing (TVD) in the scalar case 

(1.21b) TV(v n+l )<TV( V n )-, 

this trivially implies (1.19) with C — 1. 

In proving relation (1.9c) for higher order reconstructions we have used the 
assumption that for h sufficiently small there are at least r + 1 points of smoothness 
between discontinuities. Consequently we cannot apply this result to the numeri- 
cal solution v n . Nevertheless, based on heuristic analysis and extensive numerical 
experimentation, we conjecture that in the scalar case 

(1.22) TV{v n+l ) < TV{v n ) + 0(h p+l ) 

for some p > 0. 
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2. Review and overview. In (15), the first paper of this series, we present 
a second-order accurate scheme which is strictly non-oscillatory in the scalar case 

(m = 1) i.e., 


( 2 . 1 ) 


JV ( oK +1 )<W 


where N 0 (v) denotes the number of local extreme in v. This scheme is a modification 
of the “second-order accurate” MUSCL scheme [22], [4], which is total-variation- 
diminishing (TVD) in the scalar case, i.e. 


( 2 . 2 ) 


TV < TV(»"). 


In order to enforce (2.2), the slope sy (1.12a) in the MUSCL scheme is subjected 
to a so called “limiter”. Due to the operation of this limiter, the coefficient m the 
O(k) term in the Taylor expansion (1.12b) becomes discontinuous at local extrema: 
Consequently e(x) in (1.9a) fails to be Lipschitz continuous at such points, which 
leads to a loss of accuracy at local extrema. In (15) this difficulty is circumvented 
by using a modified slope sy in (1.12a) which satisfies 


(2.3) 


Sj = t/; x (xy) + 0(h 2 ), 


thus leading to a globally smooth e(x) in (1.9a). 

Although the end result is a simple technical modification of the formula for 
the slope sy, the design of the scheme in [15] invokes major conceptual changes. 
Realizing that TVD schemes, independent of their particular form, are necessarily 
only first-order accurate at local extrema, we seek a weaker notion of control over 
possible growth of total variation of the numerical solution. For this purpose we 
introduce the notion (2.1) of non-oscillatory schemes, which satisfy in the scalar 
case for piecewise smooth w 


(2.4) 


TV{E h {r) tD) < TV(w) + 0{h 2 ) 
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rather than (2.2). In [16], the second paper in this series, we show that even the 
notion of (strictly) non-oscillatory schemes (2.1) is too restrictive in the sense that 
it limits the order of accuracy to 2. To enable the design of higher-order accurate 
schemes we then introduce the notion of essentially non-oscillatory schemes (1.13), 
which excludes a Gibbs-like phenomenon but allows for the production of spurious 
oscillations on the level of the truncation error. 

Another conceptual change is the removal of the “monotonicity limiters” which 
are an essential part of TVD schemes ([30]) and may cause reduction of order of 
accuracy at some points. Our new schemes are of uniform order of accuracy r. 
The control over possible growth of the total variation of the numerical solution is 
obtained by an adaptive stencil that at each point attempts to use the smoothest 
information available. This adaptive selection of stencil is introduced to the algo- 
rithm through the reconstruction step (1.14a). The number of points in the stencil, 
independent of its orientation, is always (r+ 1). 

In [16], the second paper in the series, we investigate the stability of our new 
schemes in the scalar constant coefficient case 

u t ■+■ flu* =0, a = constant. 

The exact evolution operator (1.2) in this case is just a translation with the constant 
speed a. Therefore our schemes (1.14) take the particularly simple form 

(2.5b) t/* -1 " 1 = ~ aT ^ n )' 

Due to the adaptive selection of stencil in the reconstruction step, the scheme ( 1.23b) 
is highly nonlinear; consequently the use of the standard linear stability analysis is 
inappropriate. We demonstrate this point in [16| by choosing initial data for which 
the reconstruction algorithm selects a stencil that is biased in the “down-wind” 
direction (i.e. in the direction opposite to that of the wind); a constant choice of 
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r ) 


br i 



such a stencil is notoriously unstable. Such an instability usually exhibits itself by 
the production of increasing oscillations which starts at the highest derivative and 
propagates to the function itself. The numerical experiment in [16] shows that once 
these oscillations begin to appear on the level of the highest derivative, the adaptive 
selection of stencil in (2.5b) reacts by changing the orientation of the stencil and 
thus avoids the build up of instability. 

In [16] we also investigate the initial-boundary value problem (IBVP) for (2.5a). 
Unlike the treatment of boundaries in standard finite-difference schemes we do not 
use “numerical boundary conditions’’. Instead we modify the scheme (2.5b) by 
restricting the selection of the stencil to available information. As a result the 
scheme is biased “against the wind” at one of the two boundaries. Nevertheless, 
numerical experiments show the scheme to be strongly stable. 

In the present paper, the third in the series, we turn to consider the general 
nonlinear case. The abstract form of our schemes, (1.10) and (1.14), calls for the 
evaluation of the exact solution in the small (i.e. for 0 < t < r, r small) of the IV P 
(1.1) with the initial data R(x,v n ); the latter is a piecewise polynomial function of 
x with possible discontinuities at {iy + i}. 

When R(x; v") is the piecewise-constant function (1.11) (i.e. Godunov’s scheme), 
we can express this solution in term3 of local solutions to the Rienrann problem 


r v" z < o 

(2.7) u t + /(«),= 0, u (*,°) = | v » +i x>u - 

Wh- n R(x,v n ) is a piecewise polynomial function of higher degree we cannot in 
general express the solution of the IVP(l.l) in a simple closed form. Nevertheless, 
(see [1], [6]) we can obtain a local Taylor expansion of the solution to any desired 
order of accuracy. 

We note, however, that the step of “solution -in-the-smaU” (1.14b) is followed 
by the step of “cell- averaging” (1.14c). consequently many of the fine details of the 


J 


14 


exact solution, which may be very costly to compute, are later ignored in evalu- 
ating Vj + 1 by averaging the exact solution over (xy_i, i ) . To economize on 
the cost of our schemes it makes sense to use simplified approximate “solvers that 
carry only this information which determines the value of the cell-average, namely 
the one needed to compute a numerical flux satisfying (1.17c). The study of such 
approximate solvers is the main issue of the present paper. In Section 4 we con- 
sider the scalar case; in Section 5 we extend the scheme to hyperbolic systems of 
conservation laws. 

When we consider the reconstruction (1.9) in the context of approximation of 
functions, the assumption that u;(x) is piecewise smooth with a finite number of 
discontinuities implies that for h sufficiently small these are at least (r + 1) points 
of smoothness separating discontinuities on the computational grid. Therefore at 
any point of smoothness it is possible to select a stencil from the smooth part of 
the function. Although the x-behaviour of weak solutions of (1.1) is generically of 
this type, their time dependence allows for collision of discontinuities, as well as 
their collision with a boundary, e.g. solid walls. For points in a region between two 
discontinuities that are about to collide, no matter how small is h, there must come 
a time when there are not enough points to select a stencil of (r t- 1) points from 
the region of smoothness. Consequently a component-wise extension of the scalar 
reconstruction algorithm in [16] to vector functions may produce large spurious 

oscillations during this brief encounter. 

The elimination of such spurious oscillations has been a major consideration 
in designing the extension of our scalar schemes to systems of conservation laws. 
In Section 5 we show that this can be accomplished to a great extent by extend- 
ing the scalar reconstruction algorithms to systems vir the use of locally defined 



characteristic variables. 
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In Section 6 we describe in detail the algorithm for the solution of the Euler 
equations of gas dynamics. In Section 7 we present some numerical experiments. 

In future papers we shall present the extension of these schemes to 
two-dimensional problems and study the dependence of the computational efficiency 
on the order of accuracy of the scheme. 


/ 

J 
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3. Reconstruction. In this section we present a brief description of the 
reconstruction R(x; w) to be used in (1.14a); we refer the reader to [16], [11], [17] 
for mere details and analysis. For this purpose we introduce a piecewise 

polynomial function of x that interpolates w at the points {xy}, i.e. 


(3.1a) 


(3.1b) 


H m {xj]w) = tv(xj), 


Hm{ x;v>) = q m j+L( x ’, w ) for - 1 - X J + 1 ' 


where q m i is a polynomial in x of degree m. 

We take q m J + i to be the (unique) m-th degree polynomial that interpolates 
u>(i) at the (m + 1) successive points {x*}, *m(;) < * ^ *m(j) + m, that include 
x j and xy+i, i.e. 


(3.2a) 


9m. j+l( 2C *; tu ) = U, ( X «) ^° r - * - *"»0) + m ‘ 


(3.2b) 1 - m < i m {j) - J < o. 

Clearly there are exactly m such polynomials corresponding to the m different 
choices of t m (» subject to (3.2b). This freedom is used to assign to (xy,xy +1 ) a 
stencil of (m + 1) points (3.2) so that w(x) is “smoothest” in (x,- m (», x, m(j)+m ) in 
some asymptotic sense. 

The information about smoothness of w(x) is extracted from a table of divided 
differences of w. The latter can be defined recursively by 


(3.3a) 


U>[x,] = tu(x.) 



(3.3b) 


tu[x,, X,+fc] = (u/[x,- + i, ...,X I + fc] - w[x u .... X,>.fc_i])/(x, + fc X,). 


17 


It is well known that if w is C «, in [x,, x, +fc j then 


(3.3c) 


l a 

W \ X *I — »*»+fcJ = jfcf r « ^ £«./t < x t >fc. 


However if w has a jump discontinuity in the p-th derivative in this interval, 
0 < p < k, then 


(3.3d) u>[x,-, ..., z, +fc ] = 0(fc- fc+ P[t(;<P>]); 

here [tu (pl ] denotes the jump in the p-th derivative. (3.3c)-(3.3d) show that 
|w[x,, ...,x, +fc j| provides an asymptotic measure of the smoothness of w in (x„ x l+fc ), 
in the sense that if w is smooth in (x,,,x $1 +fc) but is discontinuous in {xi J} Xi 3+k ), 
then for h sufficiently small |u;[x lt , x, i+fc ]| < |u;[x la ,...,x lj+fc ]|. Hence the prob- 
lem of choosing a stencil of points for which w is “smoothest” is basically the same 
as that of finding an interval in which w has the “smallest divided differences.” (see 
[16], [11] for more details). 

In [11] we propose the following recursive algorithm to evaluate t m (j). We 
start by setting 

( 3 ‘ 4a ) *i (;)=;, 

i-e. + 1 is the first-degree polynomial interpolating u at x, and x )+l . Let us 
assume that we have already defined i k (j), i.e. q k} + 1 is the k - th degree polynomial 
interpolating w at 

X *kU)t •••) *•*(»+*• 

We consider now as candidates for q k+l , ] + i the two ( k + l)-th degree polynomials 
obtained by adding to the above stencil the neighboring point to the left or the 
one to the right; this corresponds to setting i fc+ ,(;) = »*(;) - 1 or t *+,(;) = t k {j), 
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respectively. We choose the one that gives a (*+ l)-th order divided difference that 
is smaller ia absolute value, i.e. 

(3.4b) 

*Wi(;) = { ” 1 lf < ! w b. fc (y),.-.,x, fc (» + fc +1 ]j 

l l k[j) otherwise. 


In [16] we analyse this interpolation technique for a piecewise smooth function 
w and show that: (i) wherever w{x) is smooth 

(3,5a) d^ Hm ^ X] = £k w ( x ) + 0(h m+l ~ k ), o <k<m; 


(>0 H m {x\ w) is an essentially non-oscillatory-interpolation of w in the 


sense that 


(3.5b) 


TV{ff m (.;w)) < TV(w) + 0(A ,n+1 ). 


We turn now to describe two different techniques to solve the reconstruction 

problem (1.9) in terms of interpolation. (See appendix for an algorithmic descrip- 
tion). 

(1) Reconstruction via Primitive function: Given cell averages w } of a piecewise 


smooth function w 


1 /*»+* 


we can immediately evaluate the point-values ol the primitive function 1V(: 

(3-7n) W-(x) = /’ w(y)iy 

* Irt 


(^y + l) — h } W } . 


w i*) = ~£ W i x ) 




we apply interpolation to the point values (3.7b) of the primitive function W{x) 
(3.7a) and then obtain an approximation to u(x) by defining 

(3.8) R{x;u>) = -^H t {x]W). 

We note that this procedure does not require uniformity of the mesh. 

The primitive function W(x) is by one derivative smoother than u>(x), therefore 
it follows from (3.5a) that wherever W(x) is smooth 

£*(*•,») 

thus we get from the definition (3.8) that 

(3.9) ^ fw(* ) + 0(h r ~ l ), 

which implies (1.9a) for / = 0. 

The conservation property of the reconstruction (1.9b) follows immediately 
from the definition (3.8): 

(3.10) [ 1+3 R{x:w)dx = ^-[H r {Xj+i\W) - H r {Xj_^ t W)\ 
h i Jx t x 

f 

[ The non- oscillatory nature of the reconstruction (1.9c) follows primarily from 

the non-oscillatory nature of the interpolation (3.5b), see [16]. 

We denote the reconstruction via primitive function (3.8) by RP. 

(2) Reconstruction via deconvolution: We assume that the mesh is uniform and 
consider the given cell- averages tZ>y to be point values of u)(z), the globally defined 
sliding-average function (1.3) of w , i.e. 

(3.11a) tD,- = fl>(x,), 
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where 


i f k/2 

(3-llb) U-kl-H,***'*' 

Expanding w{x + y) in (3.11b) around y = 0, we get 

P-m) £“ fckV ‘ ,li 


where 


(3.126) 


ctk = 


0 k odd. 

2~ k /{k +■ 1)! k even 


Multiplying both sides of (3.12a) by hfdf/iJ and then truncating the expansion in 


the RHS at 0{h r ), we get 


(3.13a) fc'tS< , '(x)= Y, e k h‘ + '«> ( * + "(^ + O(V) 

v ’ *=0 

, w <1 10-1 f rtr / — 0 r — 1 in a matrix form, we obtain 
Writing the relations (3.13aJ tor l — u, — » r 



Let us denote the coefficient matrix in the RHS ol (3.13b) by C. This matrix 
is upper triangular and diagonally dominant. Multiplying both sides of (3.13b) by 


c 1 from the left we get 



13c) 


hw • (x) 

-1 

hv'(x) 

• 

• 

- C 

• 

• 

• 

,r-l (r-1), J 


\ .1-' -tr-l 


+ o(h r ) 


W 


w 


Given Wj we interpolate tZ)(x) by H m {x\w) with m > r - 1. Since u/(x) is 
smoother than u>(x) it follows from (3.5a) that 

= £>*(*) 

wherever tu(x) is smooth. We note that although H m is only continuous at x } , the 
one-sided derivatives at x 3 ± 0 do satisfy the above relations, i.e. 

(3.14) ±0;d} ) = + C>(^ rn + 1_fc )- 

Next we define 


(3.15a) D 0 j = Wj 

(3.15b) Di,j = h l i W{~H, n {Xj - 0; tZ>), ~ H m {x } +0;ti>)) 

for 1 < / < r - 1, 


where A/(r, y) is the min mod function 


(3.16) 



min(|x|, |y|) if sgn(x) 
otherwise 


sgn(.v) -- s 


Clearly 


(3.17a) 


D,., --h'w m {z,)+0(h r )-, 
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using Dj = {Dqj, ..., D r -ij) T to approximate the vector, on the RHS of (3.13c) 
we get that 

(3.17b) Dj = C~ l Dj 

satisfies 

(3.17c) Dj = (u^x,), Au/(x J l,...,/i r-1 u/ r “ 1 *(xy)) r 4- 0(h r ). 

Finally we defu e 

r_1 1 

(3.18) £(x, tZ>) = ^2 —Dk,j[{x - Xj)/h\ k for \x - Xj\ < hj 2. 

fc= o 

We note that since C is upper triangular D/,y in (3.17b) can be computed by 
back-substitution, i.e. we set 

(3.19a) D r -i'j = D r -ij 

and then compute for k = r - 2,...,0 

r— 1 

(3- 19b) Dk. } = At.y - 'Yl a ‘ Di j- 

l=k + 1 

It follows immediately from the definition (3.18) and the relations (3.17a) that 
wherever ty(x) is smooth 

(3.20) ~R(x-w) = -^w{x)+-0{h r ~ l ); 

this for / = 0 implies (1.9a). The conservation property of the reconstruction (1.9b) 
follows from 

( 3 2»i/_ 

= D 0 .j = Ulj . 


n . x , V-i Dk , i 1 

+ «)(<» = 


r h/7 r ~‘ 

7 / = Do :] + at Dk.j 

J-hl2 . 
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The last two equalities in (3.21) follow from (3.19b) with k = 0 and (3.15a). 

The non-oscillatory nature of the reconstruction (1.9c) follows primarily from 
the non-oscillatory nature of the interpolation tf m (x;ti>); see [16] for more details. 

We note that u>(aj is the convolution of w(x) with ^(x), the characteristic 
function of a cell, i.e. 

(3.22a) tl/(x) = (ty * ^)(x) 


(3.226) 



l/h for |z| < h/ 2. 
0 for |x| > hj 2 


Hence (3.13c) is actually a deconvolution to 0(h r ). Therefore we refer to (3.18) as 
reconstruction via deconvolution and denote it by RD. 


Remark (3.1): We note that for RP with m = r and RD with m = r - 1 the 
coefficient e(x) of h r in the reconstruction error (1.9a) is discontinuous at points 
where there is a change of orientation in the stencil of the associated interpolation; 
this may happen at critical points of the function and its derivatives. Hence the 
resulting schemes (1.14) are “r-th order accurate” (see Remark (1.4)). On the other 
hand RD with m = r yield e(x) which is globally Lipschitz •'ontinuous. thus resulting 
in schemes that are r-th order accurate in a pointwise sense. This follows from the 
fact that (3.17a) is upgraded to 

f3 - 23 ) D Lj = h l w {l) (z) + 0(h r+l ) 


which has the effect of pushing the non-smoothness due to change of stencil orien 
tation in the associated interpolation to the 0(/i r+1 ) level. 
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Remark (3.2): We note that both RD with r = 2, m = 1 and RP with r = 2 are 
piecewise linear reconstructions of the form (1.12). The slope Sj for RD is identical 
to that of the “2nd-order accurate” TVD scheme in [5], The slope for RP is the 
same as that of RD except at local extrema; where s 3 = 0 for RD while for RP 


(3.24) 


w \ x i,*j+ i] if k(«y,iy + i]| < \w\xj_i , z y )| 

w [ x j-ij x ]} otherwise. 


Although RP does not “chop” local extrema as RD, the lack of smoothness in 
(3.24) results in the same loss of accuracy at local extrema. 

We note that RD with m = r = 2 is essentially the came reconstruction that 
gives the non-oscillatory second order accurate scheme of [151. 
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4. Scalar conservation laws. The abstract form of our scheme calls in 
(1.14b) for the evaluation of the exact solution in the small of the IVP (1.1) with 
the initial data R(x\v n ). This step is followed by the cell-averaging operation in 
(1.14c) which results in the conservation form (1.16). Thus we are spared the task 
of having to compute a global solution. All we need to do is evaluate 

( 4>1 ) = + »?))(*?■ 

To simplify our notation let us denote v h {x,t n + t ) by v(x, f). Thus v(x,t) is 
the solution of 


(4.2a) v t + f{v) x = 0 

with the piecewise-polynomial initial data 

r— 1 

(4.2b) u(z.O) = R( x\v n ) = ^2bjj(x- x 3 ) 'll for x } _x < x < x J+ i 

/= o 3 2 

in the time strip -oo < x < oo, 0 < t < r, r small. 

The solution v(x, t), for sufficiently small r, is composed of sections of smooth- 
ness seperated by “fans” that emerge from the discontinuities at {x f+ i}. We use 
here the term “fan” loosely, allowing a “fan” with zero spread which is just a curve. 
In the linear case discontinuities propagate along characteristic curves; in this case 
all the “fans” are just curves. In the nonlinear case the “fans” with zero spread are 
shock curves, while “fans” with positive spread are rarefaction fans - or possibly a 
succession of rarefaction fans seperated by contact shocks in the case of non-convex 
flux. We denote by v : (x, i) the section of smoothness of v(x, t) that is connected to 
the polynomial data in (xy_i, x J + i). 

A global description of v(x, <) can be quite complicated. Fortunately all we need 
is «(*y+i,0 for small t , which can be easily described in terms of v y (x, t), v j4 .i(x,t) 
and the “fan” eminating from x = x J+ i at t = 0 as follow-- If for t > 0 the “fan” 



r 

l 

\ 

i 

f 

\ stays to the right of x = x } +± then v(x 7+ ^,t) = vy(x J+ i.,t); if this “fan” stays to 

t 

\ the left of x = Xy + i then v(x ]+ ^,t) = ty +1 (xy + i,i); if the “fan” covers x = x ; + i , 

[ then u(x J + ^,t) = constant = VX0;t/y(xy + ^,0),t;y+£(xy + ^,0)). Here V(x/t-,u L , u R ) 

’ denotes the self-similar solution of the Riemann problem 

f 
r 

L f ur Z < 0 

! (4.3) u t + /(ti) x = 0, u(x,0) = < 

! {u R z>0, 

% 

i with constant U£, and Mr. We note that the “fan” covers x = x J+ x only when it contains 

a sonic centered rarefaction wave (i.e. one that includes a point for which f = 0); this 
wave retains its self-similar form as long as it does not interact with shocks. Therefore if 
we choose r sufficiently small so that no shock crosses x = x J + i for 0 < t < r, we can 
| express /(u(xy + ^,f)) by 

(4.4) 

{ /(vy(x J+ £,t)) “fan” stays to the right of x = x ;+ i 

f R ( V J (x y+ ^ , o) , 1 /y +1 (xy . 0) ) “fan” covers x = x y+ i . 

/ (uy-f i (xy + t , f ) ) “fan” stays to the left of x = x J + i 

Here f R denotes the flux at x = 0 of the solution to the Riemann problem (4.3), 

1 

\ 

\ l.e. 

(4.5) / R (ui,u 2 ) = /(V(0;ui,u 2 )); 

using the foviuuia in [23] it can be expressed by 

if '-i < u 2 

max . 1 >«>»,/(“) »f«i> u 2- 

When f(u) is a convex f’mction of u, i.e. f"(u) > 0, /(u) may have only a single 
local extremum which is a minimum; let us denote its location by u # . Using this 
fact in (4.6a) we can express / fi (uj,tt 2 ) in the convex case by 
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(4.6b) 


Here 

(4.7) 


/*(« i»tt 3 ) = 


' /(« l) 

if 

u. 

< 

Ui 

3 

VI 

/(«.) 

if 

Uj 

< 

U, 

c* 

VI 

fM 

if 

«x 

< 

u 2 

< Uj 

/(» i) 

if 


> 

u 2 

and a(ui,u 2 ) > 0 

. /(«2) 

if 


> 

U 2 

and a(u 1 ,u 2 ) < 0. 


a(u l ,u 2 ) = [/(u 2 ) - /(tii)]/(u 2 - Uj) 


is the speed of the shock with u L = uj and u R = u 2 in (4.3). We remark that 
(4.4) is deliberately formulated in terms of /(v(r jV i,<)) rather than v(z in 
order to remove ambiguity in the definition when v is discontinuous at x , i. The 

continuity of f(v) in this case follows from the Rankine-Hugoniot relation for a 
stationary shock. 

We turn now to derive a simple but adequate approximation to the numerical 
flux (4.1), which is 

^• 8 ) f } + j = ~f Q /(«(x J+ i,f))dt 

with the integrand given by (4.4). Note that the integrand is a smooth function of 

t. 

The first step is to discretize the integral in (4.8) by using a numerical quadra- 
ture 

\ f T K 

(4-9) - / g(t)dt = £ot* g{0 k r) + 0(r r ); 

Tj ° JTo 

thus 

K 

( 410 ) 

fc-f) 


L ... 
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The second step is to approximate v,(z,t] in (4.4) by its Taylor expansion 
which is obtained by the following local Cauchy Kowalewski procedure. We start 
by expressing d l v(x,0)/dx l at x y by 


(4.11a) 


dx l 


-{.* 


for 0 < / < r - 1 
for l > r. 


Next we use (4.11a) to evaluate 


(4.11b) 


d l v 


<9z* °) for all / and 0 < k < l 


by taking derivatives of the partial differential equation (1.1a) in the following 
ordered way 


(4.11c) 


r v t = -f'v x 
v « =-[/>*) 3 + /' v „] 

vtt ~[f"v t v x + f’v xt ] 

< v xxt = -[/"'( V X ) 3 + 3f"v x v xx + /'«,„] 

Wxtt = -[/'"(v*) 2 v t + /"(2v x v xt + v t v xx ) + /' V xxt ] 

Vttt = -\r{'H?v x + /"( 2v t v xt 4- v x vtt) + /'v x „] 

. etc. 


and then compute (4.11b) by successively evaluating the RHS of (4.11c); note that 
this procedure always uses known values which are either initially given by (4.11a) 
or previously computed in the algorithm (4.11c). We observe taat 

(4.12a) CJx, t) x,, 0 ) (z-z;) fe t l ~ k 

h^ dxkdtl ~ k H (/-*)! 

satisfies 


(4.12b) 
and that 


v ; (r,0) = Vj(x,0) = v(x.O) 


for r i < 
J 2 


T < X 


>M' 


U,(x,0 = v,(x,0 + 0(/T), 


(4.12c) 
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wherever Vj(x,t) is well defined. 

The last step in our derivation of the numerical flux is to approximate 
f{v{x J+ x,t )) in (4.4) by 

( 4 - 13 ) /(v(r y+ i,f))» f P {vi(x j+ ^,t) t »y+i(* y+ x,0). 

where f R is (4.5) - (4.6). The resulting numerical scheme is 
(4.14a) = -/,_*) 

/c 

(4.14b) 4+1 = 

fc=o 

In the following we show that the numerical flux / J+ i in (4.14b) is an adequate 
approximation to the “abstract numerical flux” (4.1). 

We start by proving that the scheme (4.14) is r-th order accurate in the sense 
of (1.8). To do so we take in (4.14) v" = u(xy,t n ) where u{x,t) is a smooth (either 
globally or locally) solution of (1.1) and show that 

( 4J5 ) /y+A = ; f /(u(x ;+ i,t n + rj))dr] + 0{h r ). 

T j o 

When we apply the reconstruction R to €. n we get from (3.9) and (3.20) that 
d k d k 

( 4 -l 6 a) dx* R ^ Z ' = tn ^ + °( /ir_fc ) for 0 < k < r - l. 

Consequently it follows from the Cauchy- Kowalewski procedure (4.11) - (4.12) and 
r — 0(h ) that 

(4.16b) 5,-(xy + i, 0 = u(x J + x,t n 4- t) + 0{h r ) for i = 1 

/*(ui,u 2 ) is Lipschitz-continuous with respect to and u 2 , and it is consistent 
with /(u) in the sense that f R (u,u) = /(u); therefore 

; (4.16c) /*(»!. «a) = /(«) + 0(|tt - Mil + ju - u a |). 
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Applying (4.16c) to (4.16b) we get that 


(4.16d) 


/*( w j( z /+a. 0> vy +1 (z J+ x,*)) = /(u(z J+ x,f„ + f)) + 0(h r ). 


Finally using the assumed smoothness of u(z,f) and the order of accuracy of the 
numerical quadrature (4.9) we obtain (4.15). 

Next we consider the constant coefficient case 


(4.17a) 


v t -t- av x =0, a = constant. 


Here the “fan” in (4.4) is the characteristic line 


(4.17b) 


I y+^(0 — x j+i + 


(4. 17c) vjiz, t) = v(x - at, 0) = R(x - at ; v n ) for z y _ x (t) < z < x J+ j ( t ). 

Since Vj(x,t) in (4.17c) is a polynomial of degree r - 1 in (z - at) we get that 


(4. 17d) 




this implies in (4.11) - (4.12) that 


(4.18a) 


Vj{ X ,t) = Vj(x,t). 


Hence 


(4.18b) 


«>+i(*> +4 ,0) s /(«/(*,>*, 0). 


Since the numerical quadrature (4.9) is exact for polynomials of degree r - 1 we get 
that the numerical flux / >+ i (4.14b) is identical to (4.8). It follows then that the 
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numerical scheme (4.14) in the constant coefficient case is the “abstract scheme” 
(1.14), (1.16), i.e. 

(4.18c) v” +l = R{ Xj - ar ; v n ). 


We observe that since the “fans” in the solution v(x, t) in the constant coeffi- 
cient case have zero spread, the evaluation of /(v(z J+ i,£)) in (4.4) involves only the 
smooth parts of the solution u,(z,f). The “fans” in the numerical approximation 
mark the domain of validity of the Taylor expansions vy(z,t). Therefore the only 
role of the Riemann solver in the formulation of the numerical flux (4.14b) 


(4. 19) f R (t/y [x j+ iJ). vy + 1 (z i+ 1,0) 


/(«*(*>+ *>0) if ° > 0 

/(«>+ i(*y+*.0) ifa <° 


is to serve as a pointer, i.e. to identify whether x = z J + 1 falls into the domain of 
validity of Vj or into that of vy+j. Since u,(z J+ i,0 = t»(z J+ i — a£,0), the use of 
the Cauchy-Kowalewski procedure is equivalent to that of a characteristic method 
that traces the characteristic curve through (z J+ i,0 to the initial data. 

Next we consider the scalar IVP (1.1) with convex f(u) and smooth initial data 
u 0 (z) and we show that the above interpretation of the numerical approximations 
applies to this nonlinear case as well. The numerical solution t>y ss u(z ; ,t n ) typi- 
cally forms a monotone transition of 1 - 2 points across shocks and stays close to 
u(xj,t n ) in the smooth parts of the solution, (see the numerical experiments with 
u ( uu z = 0 and u(z,0) = sin wx in section 7). Let us now examine the discon- 
tinuities of R{x-,v n ) at {r J + i} and the nature of the “fans” eminating from these 
points. Relation (4.16b) with t = 0 shows that the jump at z J + i in the smooth part 
of the solution is of the order of the local error, say 0(h p ) with 0 < p < r. Hence 
the “fan” emerging from z J+ i in a region of smoothness is either a shock curve or 
a rarefaction fan with 0(h p ) spread. On the other hand in the vicinity of shocks of 
u(z, t „ ) the si2e of this jump is 0(1); however the “fan” is necessarily a shock curve. 
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We see therefore that the global picture is very similar to that of the constant co- 
efficient case, i.e. the “fans” separating {vj(x, <)} are either shock curves with zero 
spread or rarefaction fans with 0(h P) spread (these can be thought of as “blurred” 
characteristic curves); these “fans” are the boundaries of the domains of validity of 
the Taylor expansions v,(i,f). We note that the value given by differs 

by 0{T r ) from that obtained by solving the nonlinear characterise relation for v 

(^•20) v = R(x J+ 1 — a(v)t; v n ). 

Hence the use of the local Cauchy- Kowalewski procedure is again computationally 
equivalent to tracing the characteristic curve through to the initial data. 

Since the evaluation of f(v(x J+ ±,t)) in (4.4) essentially involves only ij(x J+i .,t) 
and t >y+i(*/ + i,0« the role of the Riemann solver in the numerical flux (4.14b) 
is again that of a pointer, i.e. to identify to which domain of validity x = x i 
belongs to. This suggests that f R in (4.14b) can be adequately replaced by the 
simpler expression f ROE which corresponds to Roe’s approximate solution of the 
Riemann problem (see [25], [14]): 

f ROE {u 1 ,U 2 ) = ~[/(u,) + /(u 2 ) - [S(ui, U 2 )|(u 3 - Ui)j 

(4 21 ) _ f /(“») if a(«i,«a) > 0 

l /(« 2 ) if a(u, , u 2 ) < 0, 

where a(ui,ti 2 ) is defined by (4.7). Observe that f ROE in (4.21) satisfies (4.16c) 
and therefore the modified scheme remains r-th order accurate. 

The heuristic analysis presented above is applicable only when all the discon- 
tinuities in the solution to the IVP (1.1) are shocks; discontinuities that are not 
shocks may be present in the solution either by being introduced through the ini- 
tial data u 0 (z) or as a result of a shock-snock interaction in the non-convex case. 
Clearly f ROE in its form (4.21a) should not be used when the solution contains 
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• •*. anv discontinuity with a(u L ,u R ) - 0 as 

a sonic rarefaction wave since it admits any ai 

,, u well known and there are many ways to 
a Stationary solution. This problem >s well known 

overcome it (see (131, (26|, (9( end Section 7). 

In section 7 we present numerical experiments testing the performance 

hele ,4 ,4) in the solution o< the Riemann .VP (4.3, where /(.) U oon-convex 
scheme ( . ) . n M in others not reported here, 

and u(un , .») = 0. In all these expenments, as we 1 
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5. Systems of conservation laws. In this section we extend the reconstruc- 
tion algorithm of Section 3 and the solution-in-the-small procedure of Section 4 to 
the case of hyperbolic systems of conservation laws. 

As always we are interested only in “computable” solutions and therefore as- 
sume that the initial data u 0 (z) in (1.1b) are such that u(z,<), which is a vector 
function of m components u = («!,..., u m ) T , is a piecewise smooth function of z 
with a finite number of discontinuities. Given cell-averages u” = u(z_,, £ n ), it seems 
natural from the point of view of approximation theory to reconstruct u(x,t n ) by 
applying the scalar reconstruction R to each of the scalar component ii£, i.e. 

(5.1) R(z; u") = (R{x‘, a" ), .-, £(z; Of; 

here R denotes vector-reconstruction. However, componentwise reconstruction 
seems natural only if we disregard the time-dependence of u(x,t) which allows 
discontinuities in the solution to collide with each other. 

We recall that the scalar reconstruction is non oscillatory only if discontinuities 
are seperated by at least r-f 1 points of smoothness, where r is the order of accuracy. 
Consequently the component-by-component reconst ruct ion (5.1) may cease to be 
non-oscillatory around the discrete set of points ( x c ,t c ) where discontinuities of 
t i{x,t) interact. In the following we describe an algorithm to reconstruct u(x,t n ) 
from u n which avoids this difficulty by decomposing u n into m locally defined scalar 
characteristic variables. 

We start by examining the constant coefficient case /(u) = Au, where A is a 
constant m x m matrix 

(5.2a) u t 4- Au x = 0 


(5.2b) 


u(z,0) = u 0 (z). 
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We note that the eigenvalues {a fc } as well as the eigenvectors {r fe }, {l k } are 
also constant. We assume that 


(5.3a) 


fli < 02 < • • • < a r 


( 5 -3b) l t r } = Sij. 

We define the A:-th characteristic variable w k by 
(5.4a) iu k = l k u. 

It follows then from (5.3b) that 

m 

(5.4b) u — ^ w k r k . 

fc=x 

Multiplying (5.2) from the left by l k we see that w k {x, t) satisfies the following 
scalar IVP 


(5.5a) 

(5.5b) 

the solution to which is 
(5.5c) 


(to fc ) t 4 - a k (w k ) x = 0 
w fc (i,0) = l k u 0 ( 1 ) = w k (x), 


w (x, t) — Wq (x — a k t). 


Using (5.4b) and (5.5c) we can express the solution u(x, t) of the constant coefficient 
IVP (5.2) by 

u(x, t) = w 0 (x — a k t)r k . 

Let us now consider the following initial data in (5.2b) 



UL 

x < x L 

(5-6a) 

«o(*) = < U M 

XL < X < Xr 


{ «/? 

Xr •' X 
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First let us consider the case x L = x R = 0 which is the Riemann IVP (4.3). The 
solution u(x,t) is a self-similar solution V[x/t\ u^, up) of the following form 


(5.6b) 


tx(x,i) = V(x/t ; u l ,ur) 


Uc x ft < <Zi 

^ cifc ^ x/t ^ Ojt+i) 1 ^ ^ ^ rn 1) 

k a m < ^/£ 


where 

(5.6c) 


k 

u k — Uh 4- ^ ^ (u>p ~ , 1 < fc < m — 1. 

»=i 


In the case ip > x^, in (5.6a) the solution u(x,t), for f small, is 

(VC^itti.uw) forx<X£, + o m t 


(5.6d) 


u(x,t) 


\ U M 


for i£, + a m t < x < xp 4- a^. 


( V r (S^n. ;UM) U/ j) for xr + a.it < x 

As t increases, the discontinuity in the fc-th characteristic field originating at x = 
xl will eventually collide with any discontinuity in the /- th field, l = 1, ..., k — 1 
originating at x = xp. 

The example (5.6) demonstrates the difficulty encountered in using the com- 
ponentwise reconstruction (5.1). We may get oscillations for small t in both (5.6b) 
and (5.6d) since the discontinuities are too close due to the self-similar nature cf 
the solution to the Riemann problem. Later on we may get more oscillations in 
(5. fid) as discontinuities collide. 

We observe that there are no such problems with u/ k (x,t) — w k {x - afct). 
Therefore it makes sense to use the scalar reconstruction R{x\ u> k ) to define 

m 

(5.7a) R(x;u) = ^2 R ( x '^ k ) r k 

k = \ 


where 


(5.7b) 


= M- 


\ 
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We generalize (5.7) to the nonlinear system case by using locally defined char- 
acteristic variables. To reconstruct u from u in (xy_i, Xy + i) we define 

m 

(5.8a) R(x; u) = £ #(x; u> fc (uj))r fc (uy) for x.,_i < x < x j+ 1 , 

where the mesh function w k (Uj) = {wf(wy)} is defined by 

(5.8b) tD?(fij) = /*(%)«, for j-p<i<}+ p; 

here p is the desired order of reconstruction. 

In Section 7 we present calculations for the Euler equations of gas dynamics 
with the initial data (5.6a). The results of these calculations (as well as those of 
shocks reflecting from a wall) demonstrate that the reconstruction (5.8) works well 
also in the nonlinear case. 

We turn now to describe our scheme in the case of hyperbolic systems of con- 
servation laws. This scheme is identical in form to (4.14): 

(5.9a) v” +1 = v? - A (/ /+ 1 - /y_i) 

K 

(5.9b) / i+ 1 = ^afc/ H («i(*y+*./3fc0. 5 >+i(*j+*i/ ? fc r ))* 

fc = 0 

The derivation of (5.9), although different in some details, is basically the same as 
the one presented in Section 4 for the scalar case. Rather than repeating ourselves 
we shall use the formulae of Section 4 (which are to be interpreted here in a vector 
sense), and point out the differences whenever they do exist. 

The problem to be solved in the “solution-in-the-small” step of the algorithm 
(1.14b) is as before (4.2). The general structure of the solution v(x,f) is similar 
to that of the scalar case, i.e. it is composed of sections of smoothness seperated 
by “fans” eminating from the discontinuities at {xy + i}. As in the scalar case 
we can use a local Cauchy-Kowalewski procedure to approximate u,(x,<), the sec- 
tion of smoothness of v(r,£) that is connected to the polynomial initial data in 



i 


w 
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by Vj{x,t) in (4.12) to any desired order of accuracy. Since /(u) is 
now a vector, /'( u) is a matrix, /"(u ) is a tensor and so on; consequently (4.11c) 
has to be replaced by a much more complicated expression. Rather than doing 
this we shall present in Section 6 an algorithm to carry out the Cauchy- Kowalewski 
procedure in the specific case of the Euler equations for gas dynamics. 

Next we consider the “fan” that eminates from the discontinuity at x j+ a. As 
in the scalar case this “fan” starts at t = 0 as a self-similar solution to the Riemann 
problem (4.3), which in the system case is a packet of m fans corresponding to the 
different characteristic fields. A major difference from the scalar case is that (except 
when the initial data in (4.2b) are piecewise constant) the “fan” emerging from x J+ 1 
at t = 0 immediately loses its self-similar nature. Therefore it is no longer possible 
to express t>(xy + i,f) in a simple closed form as we did in (4.4). However v(xj±L } t) 
can be expressed to any desired order of accuracy via a local Taylor expansion of 
the various curves in the “fan” and the states in between (We refer the interested 
reader to (1] where Ben-Artzi and Falcowitz describe such an expansion for the Euler 
equations of gas dynamics). Thus as in the scalar case, although at a considerably 
more effort, it is possible to obtain an explicit expression that approximates the 
“abstract numerical flux” (4.8) to any desired order of accuracy. 

We turn now to show that the numerical scheme (5.9) is an adequate approxi- 
mation to the “abstract scheme” (1.16). First we observe that relations (4., 6) hold 
also for the system case; therefore (4.15) follows in exactly the same way as in the 
scalar case and consequently the scheme (5.9) is likewise r-th order accurate. 

Next we consider the scheme in the constant coefficient case (5.2). Since both 
the PDE (1.1a) and the scheme (5.9) decouple into m scalar relations for the char- 
acteristic variables w k in (5.4a), we can apply the analysis of the scalar constant co- 
efficient case to systems in a characteristic-wise fashion. It follows then from (4.17) 
- (1.18) that the numerical flux (5.9b) is exact and that the numerical scheme (5.9) 


39 


is identical to the “abstract scheme” (1.16). Let us examine now the structure of 
the solution v(x,t) : The “fan” eminating from iy + i has the same form as (5.6b) 
except that u, k and ur are now functions of x and t. The section of smoothness 
Xy_ i + a m t < x < x J+ i + ayt is also the domain of validity of the Taylor expansion 
ij(x,t). We note however that l k tiy(x,£), which is the Taylor expansion of w k (x,t), 
is valid in the larger domain Xy_i 4- a k t < x < Xy + i + a k t. Next let us examine 
the role of f R in formulating the numerical flux (5.9b): 

(5.10a) /*(5i(*j + i,0.«y+i(*y+A.0) = J2^ + ^ lk 

k 

Vj+i ( x j+$ > 0] r fc» 
k 

where 

(5.10b) (a fc ) + =max(0,a fc ), (a k )~ - min(0, a k ). 

We see from (5.10) that as in the scalar case the role of f R is that of a pointer, 
i.e. to identify for each characteristic variable u) k — l k v to which domain of validity 
of {l k v,} does x — x J+ i belong to. Since l k t)i(xy + i,£) = l k v(x J+ i - a k t, 0), the 
use of the Cauchy- Kowalewski procedure in this fashion is again computationally 
equivalent to that of a characteristic method. 

In the following we argue that except for the discrete set {(x c ,t c )} of interac- 
tions, the above interpretation can be applied to the nonlinear case as well. Unlike 
the scalar case we do not consider in this paper the “non-conv< x case” for systems 
and assume that each characteristic field is either genuinely nonlinear or linearly de- 
generate (see [19]). When we consider the IVP (4.2) in the context of the numerical 
scheme where v(x,0) = R(x \v n ) we see that the “fans” in the solution v(x,t) are 
related to the global structure of u(x, t n ) in the following way (see figure 14 and fig- 
ure 16): When u(x.£ n ) is smooth, the “fan” has the basic structure of the constant 
coefficient case linearized around u(iy + i,£ n ), except that the k - waves may have a 
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spread of 0(h p ). When ij+i is in the vicinity of a shock of u(z,t n ), the fan is 
essentially a shock wave with small perturbations in the other fields. We see that 
typically (excluding interactions) the “fan” eminating from z J+ 1 in the solution 
t;(x, t) is degenerate in the sense that except possibly for a single large shock (or a 
contact-discontinuity) all the waves in it are weak. This heuristic analysis suggests 
that /(v(z J + i,<)) can be adequately approximated by a local Roe’s linearization; 
this linearization is exact for a single shock or a contact-discontinuity and amounts 
to a characteristic approximation for weak waves. 

As in the scalar case, f ROE is obtained by a local linearization with respect to 
a particular average u = u(tin,Ufl) for which 


(5.11a) 


/(u«) - /(*l) = - u l)- 


fROE j s defined as the flux at x — 0 of the solution to the constant coefficient 
Riemann IVP: 

u t + A(u)u x = 0 

, n , f ul x < 0 
“(z.Oj = S . n 

l u R x>0, 

which can be expressed as 

1 m 

(5.11b) J ROE {u l ,u r )= -[/(tin) +/(u R ) - X^ fc ( Ut > u *)M“)i r fc(«.)] 


k= t 


where 


(5.11c) 


M u 4>“r) = /fc(«)(«R - «t); 


here a fc (u), l fc (u) and r fc (ti) are evaluated with respect to the Jacobian matrix A(u). 
The derivation of Roe’s Riemann solver is well documented in the literature (see 
[25], [8], [9], [14]). In Section 6 we describe f ROE for the Euler equations of gas 
dynamics 




niui.y#p iJii 
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Finally let us examine the performance of the scheme (5.9) during an inter- 
action of discontinuities in the solution u(x,t) of the IVP (1.1). We observe that 
it takes some time until the outcoming waves can be properly described on the 
computational grid. Till then R{x) v"), which is based on polynomial interpolation, 
can only be a crude approximation to u(z, t n )- Under these circumstances we ex- 
pect the “fans” in the solution v{x,t) (4.12) that originate from discontinuities in 
the interaction zone of u(x,f„), to be adequately approximated by the self-similar 
solution to the local Riemann problem. We note that once the outcoming waves 
are properly resolved on the computational grid, the previous analysis applies. 

In Section 7 we present numerical experiments where the scheme (5.9) with 
f R replaced by f ROE (5.11) is applied to an interaction problem for the Euler 
equations of polytropic gas. In all these experiments the scheme (5.9) has developed 
the correct structure of the solution. 

We remark that the scheme (5.9) with j n form (5.11b) admits a 

stationary “expansion shock” as its steady solution. This can be easily rectified 
by adding entropy viscosity terms for the genuinely nonlinear characteristic fields. 
(See [13], [14], [9]). 
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6. Euler equations of gas dynamics. In this section we describe how to 
apply the scheme (5.1) to the Euler equations of gas dynamics for a polytropic gas: 

(6.1a) ti« +/(«)«= 0 


(6.1b) 


u = (p, rn,E) T 


(6.1c) 


/(u) = qu 4- (0 ,P,qP) T 


(6. Id) P = - 1)(E - 

Here p, q, P and E are the density, velocity, pressure and total energy, respectively, 
m = pq is the momentum and q is the ratio of specific heats. 

The eigenvalues of the Jacobian matrix A(u) df /3u are 


(6.2a) 


ai(u) = <7-c, a 3 («) = <? a 3 (u) = <? + c 


where c = (^P/p) 3 is the sound speed. 

The corresponding right-eigenvectors are 


1 


1 


(6.26) 


n(u) =1 q-c I ,r 2 (u)= I q I ,r 3 («) = 
H-qc) \k<1 2 


1 

q -t c 
H 4 - qc 


here 


(6.2c) 


H = {E + P)/p = c 7 /{l - + 



is the enthalpy. 

To compute (l t (u)} which is bi-orthonormal to (r fc (u)) in (6.2b), we first form 
the matrix T(«), the columns of which are the right-eigenvectors in (6.2b) 


* ulL 


] 

J 


T(u) - (ri (u), r 3 (u), r 3 (u)) 
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and then define l k (u) to be the k - th row in T _ 1 (u), the inverse of T(u). We get 


(6. 2d) 

M u ) = 3(^2 + q/c, -b 1 q - l/c,bi) 
< l 2 (u) - (1 -6 2 ,6i9,-<>i) 

. ^( u ) = 2(^2 ~ 9 /c, -6 j 9 + l/c,6i) 

where 


( 62 e) 

cr- 

II 

T 

1 

'cT' 

to 

( 6 . 2 f) 

b 2 = ~ 9 2 6 i- 


Given {v?}, approximation to {u(rj,f n )}, we use (6.2d) - (6.2f) to evaluate 
the locally defined characteristic variables (5.8b) 


(6.3a) 


w t («? ) = l k{v")v? for i = j - r , ..., ; + r and Jk = 1, 2, 3. 


Next we apply our scalar reconstruction algorithm to each of the locally defined 
characteristic variables in (6.3a). The scalar reconstruction R{x\ u>) is described in 
an algorithmic form in an appendix; the output of this algorithm is in the form of 
the finite Taylor series in (4.2b). Thus we get for each characteristic variable in 

(6-3b) A(*;U>*(i,-))=S 6 ii(*-*>) , /« 

/ = 0 

Rearranging terms we can express the vector reconstruction (5.8a) by 

r — 1 

( 63c ) R(x;u n ) = ^bjAx - x ; )'//! 


where 

(6.3d) 




Note that wherever the solution is smooth 


b yi = JIi(P' m ’ E ) T \ x=z +0{h r ~‘) for 0 < / < r - l. 
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We turn now to describe the Cauchy- Kowalewski procedure (4.4) - (4.5) for 
the Euler equations of gas dynamics. We start by using the reconstruction output 
(6.3c) to define (4.4a), i.e. 

(6.4a) g<t, ( g >»°) = | £/./ for 0 < / < r - 1 

dz l 1 0 for / > r. 

We find it convenient to express d l v(xj,0)/dx k dt l ~ k in terms of derivatives of 
the 4-vector Z = (p,m,P,q) T . For this purpose we use (6.4a) and the relations 

m — pq 

P = (7 - 1)(£ - 

to find the x-derivatives of q and P, by 


f ~ QPx 4" PQ x —^Qx — ( m x ~ QPx)/p 
\ P x = (l ~ 1)[ E x - ±( q x m + qm x )] 



(6.4c) ( m * x ~ P<izx + 2px<Jx + W xx =* = ( m « ~ QPxx ~ 2 q x p x )/p 

\ Px x = h ~ 1 )\E XX - ~{mq xx + 2 q x m x + qm xx )] 

and so on. Having evaluated d l Z{x } , 0)/dx l for 0 < / < r - 1, we proceed to obtain 
the rest of the derivatives d l Z(xj,0)/dx k dt l - k , 0</<r-l,0<jfc</by 
differentiating the PDE’s 

(6.5a) 

(6.5b) 

(6.5c) 

and the algebraic relation 
(6.5d) 


Pt + m x = 0 
m, + ( qm) x + P x = 0 

Pt + qPx + ~iPux = o 

m=qp 



in the following ordered way: 
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Compute Z t (xj, 0) from 


(6.6a) 


' p t 4 - m, = 0 
m t 4 - q x m + qm x + P x = 0 
p t + qP x + "fPq x = 0 
. PQt + p t q = m t 


compute Z xt (xj, 0) from 


' Pxt + m xx = 0 

4 - 2 q x m x + qm xx 4- mq xx + P xx = 0 
^Xt + 9-Prx + 7- p ?xx + (1 + 7 )q x P x = 0 

, pq x t + p x qt + ptq x + qp x t = m xt 


compute Z tt (xj, 0) from 


(6.6c) 


ptt 4 - Tn x t — 0 

171 tt + qm xt 4 - mq xt 4 - P xt + ( q x m t 4- q t m x ) = 0 

D i ~ D t n _ i f . n n \ » 

* tt -ryt xt -r 11 q xt -+- \q t r x -+- 7 r t q x ) = u 
, pqtt + 2 p t q t 4 - qptt ~ m tt 


compute Z xxt (x } , 0 ) from 


( Pxxt 4 * — 0 

+ qrn X xz + mq xxx 4- P xxx 4- 3 q x m xx 4- Zq xx m x = 0 
^xxt + 9-Pzxx + iPqxxx 4 - (2 4 - l)q x Pxx + (14- 2~i)q xx P x .= 0 ’ 
' ^7**t + 2(/> x g xt 4- p xt <7 x ) 4- qp xxt 4- p xx qt 4- p t q xx - m xxt 

compute Z xtt (x } , 0) from 


' Pxtt f m xxt = 0 

m xtt + qm xxi 4- mq xxt 4 - P xxt 4- 2 (^ x fn x j 4- m x <7 x j) 

< +(m t q xx 4- ? t »TJxx) = 0 

Pxtt + + 1 iPqxxt + (7 + 1 ) (<7 x t P x 4 - q x P rt ] 4 - q t P xz 4 - 7^«<7xx = 0 

■ /*7*« + + p x qtt + q x ptt 2 {pxtqt + Ptq z t = ^»x« 


(6.6e) 


1 
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compute Z t tt(xj, 0) from 


(6.6f) 


Pttt + m x tt — 0 

m ttt + qmxtt + + Pxtt + + rn t q xt ) 

+(?2 m tt + mxQtt) = 0 

Pttt + qPxtt + 7 Pq x tt + IQxPtt + 2{q t P x t + iPtqn) + qttPx = 0 
pqttt + qpttt + 3 (ptqtt + qtPtt) = m e« 


and so on. 


We note that one can differentiate the algebraic relation (6. Id) in order to 
obtain d‘E( x 3 , 0)/dx k dt l ~ k in terms of the already evaluated derivatives of P,q and 
m, and use the derivatives of the conserved quantities p,m,E to compute Vj(x,t) 
in (4.5). However, it is more convenient to evaluate the flux f(u) and f R {u i,U 2 ) in 
terms of p,q and P\ since Vj{x,t) is smooth and the scheme (5.1) is in conservation 
form we do not really have to worry about relation (4.12b). For this reason we use 
the first, third and fourth components of Zy(x, t ) 

- % ^^a'z(z,,o) (x-z J ) fc t l ~ k 

( 6 - 7 ) z y k\ ( i-k)\ 

1 = 0 k = 0 v 

to define pj[x,t), Pj{x,t) and qy(z,f), respectively. 

Once we have computed (6.7) we can compute the numerical flux / J + i in 


(5.1b). 


An exact f R {ui , u 2 ), i.e. 


/*(« l, u 2 ) = f[V (0; u n u a )) 


where V(z/<;u 1 ,u 2 ) is the exact solution of the Riemann problem for the Euler 
equations of gas dynamics can be computed through an iterative algorithm. This 
algorithm is rather complicated, and we refer the reader to [5], [3| and [28] for its 


details. 
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To compute f R = f ROE in (5.11b) and (5.9b) all we need is to describe the 
particular average u(u!,u 2 ) for the Euler equations of gas dynamics (see [25J, [9]). 
To do so we denote the arithmetic mean of 6(u) with respect to ui and u 2 by 

(6.9a) (b) = ^(6(u 1 ) + 6(u 2 )l 

and define 

(6.9b) 4=<>/?«)/(vW. H=(SpH)/(s/p), ^ = (-7 " l) 4 

here H is the enthalpy (6.2c). Having prescribed q , H and c, we have all the 
quantities needed to define the eigenvalues and eigenvectors in (6.2). 

REMARK (6.1): The importance of using the particular average (6.9) rather than 
a simpler one is that when (u t , u 2 ) corresponds to a single shock or a single contact 
discontinuity in the solution of the Riemann problem V{x/t\ui,u 2 ), then f ROE is 
exact, i.e. 

(6.10) /*° £ («i,u 2 ) = f R {V{ 0;ui,u 2 )). 




4 
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7. Numerical experiments. In this section we present results of several 
computer experiments with the schemes (4.14), (5.9). These schemes will be referred 
to as r-th order ENO schemes (or Vth order” when applicable - see Remark (1.2)); 
ENO stands for Essentially Non-Oscillatory. 

The ENO schemes are highly nonlinear and consequently do not easily lend 
themselves to rigorous analysis. At present we have completed the analysis of the 
non-oscillatory interpolation H m (3.1) - (3.5) and have acquired a fairly good under- 
standing of the reconstruction R(x ; tZ>) ; these reconstruction results can be extended 
to a single application of the “abstract scher. e” (1.10) to piecewise smooth data. 
Unfortunately we have not been able as yet to analyse rigorously the crucial question 
of accumulation of error. Under these circumstances, computer experiments have 
become our main tool of analysis. We have performed a large number of numerical 
experiments with initial data ranging from random noise to smooth functions. We 
have studied two notions of “stability”: (i) Boundedness of a refinement sequence 
h —* 0, t — O(h) for 0 < t < T. (ii) Boundedness of the numerical solution 
as n — ► oo with fixed h and r. In all our experiments 4 we have found the ENO 
schemes to be stable under a CFL restriction of 1 and strongly so, in the sense that 
they strongly damp high frequency noise - this is probably due to the cell-averaging 
step (1.14c). 

In [15], the first paper in this series, we have presented numerical results which 
compare the second order ENO scheme based on RD with r = 2 to a “second order 
accurate” MUSCL-type scheme, which is computationally equivalent to the “second 
order” ENO scheme based on RP with r = 2. 


4 The only exception where we had to reduce the CFL number is for the initial data 
of the mesh oscillation function vj - (-1) J . This choice of initial data forces the 
ENO scheme to become linear, for vj = where 9 } is a positive random 

number, the scheme is again stable under a CFL restriction of 1: (see [ 16] for more 
details). 
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In [16], the second paper in this series, we have presented numerical experiments 
that verify our statements about the accuracy and non-oscillatory nature of the 
reconstruction R(x; w), and demonstrate the stability of the ENO schemes in the 
scalar constant coefficient case for both the pure IVP and the mixed initial-boundary 
value problem (IBVP). 

In this paper we present a sample of our numerical experiments for the non- 
linear scalar case and the Euler equations of gas dynamics in ID. The purpose 
of this presentation is to address the open questions that we could not fully an- 
swer by analysis: The accumulation of error, the adequacy of the “solution m the 
small” procedure, consistency with entropy inequalities and the effectiveness of the 
characteristic-wise reconstruction for systems. We have performed most of the nu- 
merical experiments for r = 1, 2, 3, 4, 5, 6. Since it is not practical to present six sets 
of data for each problem we usually compare r = 2, which is the current state of 
the art scheme, to r = 4 which seems to be optimal for smooth solutions. However 
presentation of a comprehensive efficiency study is deferred to a future paper. 

A. Scalar conservation laws. 

Al. Convex /(u) with smooth initial data: In this sub-section we show results 
of applying the ENO schemes to 

(7.1a) u t -t- (u 2 /2)i = 0 



(7.1b) «(*,<)) = a +/3sin(irr + Tf), 

for -1 < x < 1, ( > 0. In these calculations we have used the ENO schemes with 

f R replaced by f ROE (4.21) (without any entropy correction). 

Let Z{x,t) denote the solution of (7.1a) with 2(i,0) = sin irx, i.e. (3 = 1, a = 
7 = 0 in (7.1b). The solution Z(x,t) is smooth for 0 < t < 1/jt; when t = l/?r a 
shock develops at x ~ ±1, and stays there as a stationary shock for t > \/n. Some 
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time after its development, this shock starts interacting with the expansion wave in 
-1 < i < 1 ; this brings about a fast decay of the solution. The “exact” solution 
presented in the following is computed in 0 < x < 1 by using Newton-Raphson 
iterations to solve the characteristic relation 

(7.2) Z — sin ir(x - Zt); 

Z in (- 1 , 0 ) is obtained from Z in ( 0 , 1 ) by Z(—x, t ) = — Z(x,t). The general solution 
of (7.1) is computed from Z(x,t) in (- 1 , 1 ) by 

(7.3) u(x,t) = a + (3Z(x - at + 7 , fit). 

In Tables 1 and 2 and figure 1 we present computation of ( 7 . 1 ) with a = 1 , /? = 
1/2) 7 = 0 , i.e. u(x,t ) = 1 + %Z(x — t, |t); thus the shock develops at t — 2 /tt. 
The results are presented at t = 0.3 when the solution is still smooth. We divide 
(-1,1) into 7 equal intervals and define 

(7-4) xj = -1 + (; - l/2)h, h = 2/ J, 1 <j<J 

First we consider the pure IVP for (7.1), i.e. periodic boundary conditions at 
x = ±1. In Figures la and lb we show the results of the ENO schemes with RD at 
t = 0.3; figure la shows the second order ENO scheme, while figure lb shows the 
fourth order one. Both calculations were performed with J = 10 and CFL = 0.6. 
The continuous line in these figures is the exact solution; the circles represent the 
values of R{ij ]v N ). In Tables la and lb we list the Loo-error the L { -error at 
t = 0.3 of a refinement sequence 7 = 8, 16, 32, 64, 128 for r = 1,2, 3, 4, 5 with 
CFL = 0.6. Table la shows the results of the ENO schemes with RP while Table 
lb shows the ones with RD. The value of r c in Tables 1 and 2 is the “computational 
order of accuracy” which is calculated by assuming the error to be a constant times 
h r '\ this definition is meaningful only for h sufficiently small. 
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In figures 2a and 2b we use the same schemes as in figure 1, but with J = 16, 
and show the results at t « 2/ir (after 17 time-steps) which is the time of the 
formation of the shock. In figures 3a and 3b we show the reconstruction R(x;v n ) 
corresponding to the numerical solutions of figure 2. The squares in figures 3a and 
3b mark the values of #(xy + i ±0;v n ). R(x:v n ) is piecewise-linear in figure 3a and 
piecewise - cubic in figure 3b. 

Next we consider the IBVP for (7.1); since the characteristic speed for (7.3) 
with a = 1, /? = ^, 7 = 0 is positive, we prescribe 

(7.5a) u(-l,0 = sr(0; 

x = +1 is an outflow boundary and no condition is prescribed there. To be able 
to compare with the periodic problem we take g(t ) in (7.5a) to be the value of the 
periodic solution at x = — 1, i.e. 

(7.5b) = i + 

The point of view that we have taken in treating boundary conditions is consistent 
with the presentation of the “abstract scheme” (1.10), (1.14) as a sequence of global 
operations. Thus in the reconstruction step, as in the pure IVP case, we use the 
given cell-averages {v"}, 1 < j < J, to get R(x;v n ) for -1 < x < 1; in presence 
of boundaries we restrict the choice of stencil to available information by imposing 
the condition 

(7 6) 1 < tfc(j) < J - r for 1 < fc < r 

in the algorithm (3.4). Note that we do not use the given boundary data g(t) 
(7.5a) in the reconstruction step. The boundary data is incorporated into the 
scheme on the PDE level by considering the solution-in-the-small step to be an 
IBVP. Obviously the resulting scheme is biased “against the wind” near x — — i ; 
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nevertheless, numerical experiments in the nonlinear case as well as in the constant 

coefficient case (see |16|), indicate that the ENO schemes are stable. We observe 

that a similar choice of stencil occur, near discontinuities in the interior of the 
computational grid. 

In Tables 2a and 2b we repeat the calculations in Table 1 for the IBVP ( 7 . 1 ) 
with (7.5). In figures 4a and 4b we show the calculations of the ENO schemes based 
on RP with r = 2 and r = 4, respectively, for the IBVP (7.1) with a = 0, 0 = 
1, 7 = s. Here the boundaries x = ±1 are characteristic, and a stationary shock 
develops at x = 0 at I = l/rr. In these calculations we have treated x = as an 
inflow boundary and specified 

u(-l,f) = 0; 

r = +1 was treated as an outflow boundary. The results show the numerical 

solution with 7 = ,6 and CFL = 0.6 at t = 0.6, at which time the solution has 

already started to decay considerably due to the interaction of the shock with the 
expansion waves. 

(A2). Rtemmn IVP for non-convex /(«): In this subsection we show results 
of applying the ENO schemes to the Riemann IVP 


(7.7a) 


u t = /(«), = o, u(x, 0) 


x < 0 
x > 0 


where f(u) is the non-convex functi 


(7.7b) 


/(u) = -(u 2 - l)(u 2 -4). 


We recall that the main difficulty in justifying the approximation (4..3) is when 
m (4.4) covers x = r J+ 1 ; the same difficulty is encountered in justifying 
the use of /«>* (4.21) instead of the exact flux of the Riemann problem (4.6). 
Therefore we present two cases in which aK,u B ) = 0; in the first case x = o 
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„ covered by a centered sonic rarefaction fan, while in the second one there is a 
stationary (sonic) shock at i = 0. 

, n each case we present two sets of experiments. In the first set we use the 
ENO schemes with the exact /« which is defined by (4.6a); these results, which 
we consider to be rather pleasing, are presented in Figures 5 and 7. In the second 
set of experiments we use the ENO schemes with /“ replaced by the follow, ng 

modification of f ROE in (4.21) 


(7.8a) 


/ RO£ (u l ,uj)= -(/(“,) + /(“a) ~ max(|d(“i,“ 2 )l»c )(“2 “ill- 


The addition of the linear viscosity term -c(“ 2 - “,)/» for W < «. - the simplest 
but crudest entropy correction of (4.2.). We note that « = 0 in (7.8a) corresponds 
to (4.21), while s = 1/A (A = r/fc) corresponds to Lax’s first order scheme (18]; since 
(7.8a) satisfies relation (4.16c) the modified scheme remains r-th order accurate. In 

our calculations we take 


(7.8b) 


e = 0.1/A. 


Analysis presented in |24| shows that using (7.8a) - (7.8b) in the *second-order” 
TVD scheme of |9i results in a scheme which converges to entropy correct solutions 
for convex /(«), provided that A is sufficiently small; numerical experiments in the 
convex case |9] and the non-convex case [32| seem to verify this statement even for 
CFL number close to 1. 

The numerical results of the ENO schemes using (7.8a) ■ (7.8b) are shown in 
Figures 6 and 8. These results show that the ENO schemes converge to entropy 
correct solutions; however, the quality of the numerical approximation depends 
strongly on the formal order of accuracy of the scheme. 
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jr remar k that an entropy correction to f ROE which is more appropriate for 

| the non ‘ conv ex case is obtained by using in (7.8a) £ = e( Ul , u 2 ) which is defined by 

S' ( 7 - 8c ) £ = max[0, a(ui,ii 2 )-a Ll a R -a(u 1 ,u 2 )] j 

■1 

f 

l where j 

i 

r | 

( 78d ) a L= min a(u 1>V ), a R = max o(v,u 2 ); • 

; j 

see [13]. In this case the modified f ROE becomes computationally equivalent to the 

exact f R . 

Our purpose in presenting numerical experiments with the crude entropy cor- ! 

rection (7.8b) rather than the more appropriate one (7.8c) - (7.8d) is to demonstrate 
that the importance of the Riemann solver in the formulation of the ENO schemes 

is decreasing with increasing order of accuracy. When r = 1 R( X \ v n ) is piecewise- ] 

| constant and all the variation of the solution is contained in the discontinuities of j 

j reconstr uction. Consequently the Riemann solver is the only mechanism to j 

| describe time evolution. For r > 1, the smooth polynomial variation in the cell j 

(which is 0(h) in regions of smoothness) is generally larger than the variation in 
the discontinuities of the reconstruction (which is 0(h r ) in regions of smoothness) 

I i - see Figures 3, 14 and 16. Therefore the time evolution of the smooth polyno- 

mial part, namely the Cauchy- Kowalewski procedure, is in general more important 
than the Riemann river. The only exception is in the first few time-steps needed 
to introduce intermediate states in the solution to the Riemann IVP (7.7) where 
(ul , u R ) is not a shock. 

In all the calculations presented in this subsection we have used the ENO 
schemes with RD and CFL = 0.8. 

L. ... J 
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Case: (i): U£, = 2, uj* — 2. 


The exact solution in this case is (see figure 5a) 


(7.9a) 


u(x, t) 


'2 x/t< -0.5281529 

< g(x/t) \x\/t < 0.5281529; 
—2 x/t > 0.5281529 


here g{x/t) is a centered rarefaction wave: g{y) is the solution of 


y = f'(9) 


in the concave part of / which is 


N < n/Vg; 


^(±0. 5281529) = T0.2152504. 

In figures 5b, c, d we show the results of the ENO schemes using the exact f R 
as defined by (4.6a) for r = 1, 2, 4, respectively; in these calculations we used J = 
40 in (7.4) and N = 80 time-steps. The exact solution is shown by the continuous 
line; the circles mark the values of *(.,-;«"). We observe that the structure of the 
solution in these calculations has developed at the correct rate; this is evident from 
the fact that the location of the computed shocks is accurate. In figure 5b we notice 

the “dog-leg” which is typical of Godunov’s scheme. 

In Figures 6a, b, c we repeat the calculations in Figures 5b, c, d but with 
f n replaced by f R()E (7.8a ) • (7.8b). From these figures we see that the scheme 
develops the correct structure of the solution, but not at the correct rate. This is 
due to the fact that £ = 0.1/ A represents a fan which is much narrower than the 
initial fan in the exact solution. The location of the computed shocks lags behind 
the correct location by 8 cells for r = 1, 3 cells for r = 2 and only one cell for 
r - 4. To verify that the numerical approximations converge to the entropy correct 
solution we refine the mesh by a factor of 2 and repeat the calculations of Figures 
6a. b, c with J = 80, N = 160; the results of these calculations are shown in figures 
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6d, e, f. Since the number of cells by which the computed shocks lags behind the 
correct location remains the same, we conclude that the numerical approximations 
indeed do converge to the entropy correct solution. 

Taking into account the crudeness of the entropy correction (7.8b) we consider 
the performance of the 4-th order scheme in figures Oc and 6f to be surprisingly 
good. 

Case (ii): = -3, ur = 3. 

The exact solution in this case is (see figure 7a) 

—3 x/t < —19.5 

g{x/t) -19.5 < x/t < 0; 

—g(—x/t) 0 < x/t < 19.5 

3 19.5 < x/t 

here g(y) is the solution of 

y = m 

in the convex part of / which is |u| > \/5/6. Note that the solution (7.9b) is 
discontinuous at x = 0; g( 0) = 

In Figures 7b, c, d we show the results of the ENO schemes using the exact f R 
as defined by (4.6a) for r = 1,2, 4, respectively; in these calculations we used J = 
40 in (' 4) and N = 20 time-steps. We observe that the stationary shock at i = 0 
in these figures is perfectly resolved. 

In figures 8a, b, c we repeat the calculations in Figures 7b, c, d but with f R 
replaced by f ROE (7.8a) - (7.8b). Since the rarefaction fans in this case are not 
f sonic, the quality of the numerical approximation of the rarefaction wave in figures 

8a, b, c is similar to that of the corresponding one in figures 7b, c, d. We observe 
that the stationary shock at x = 0 in figures 8a, b, c is somewhat smeared - this 
is due to the fact that the Riemann solver corresponding to (7.8b) places a fan of 
the size \x/t\ < e around x = 0. Nevertheless, if we compare the results of the 4-th 

Li 
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order ENO schemes in the two experiments, we find that the results in figure 8c are 
only slightly inferior to those of figure 7d. 

B. Euler equations of gas dynamics. In this subsection we present numer- 
ical experiments with the ENO schemes for the Euler equations of gas dynamics 
for a polytropic gas with 7 = 1.4 (see section 6). In all these calculations we have 
used reconstruction via primitive function (RP) and f ROE (5.1l), (6.9) without any 
entropy corrections. 

(Bl). Riemann problems. In figures 9 and 10 we show the results of applying 
the ENO schemes with r = 2 and r = 4, respectively, to the Riemann problem 
(7.7a) with the initial data 

(7.10a) {pl,<1l,Pl) = (1,0,1); (p R , q R , P R ) = (0.125,0,0.10). 

In these calculations we have used the characteristic reconstruction (5.8), (6.3) with 
100 cells, h = 0.1, CFL = 0.8 and 50 time steps. 

In figure 11 we repeat the calculation of the “4-th order” ENO scheme in figure 

10 but with component-wise reconstruction (5.1). Ccmparing figure 10 with figure 

11 we see that there is some “noise” in the component- wise reconstruction which is 
eliminated by using characteristic reconstruction. We note however that the level 
of “noise” in figure 11 may be considered acceptable for practical calculations. 

The initial data (7.10a) are those of the Riemann problem proposed by Sod 
in [23], which has become a standard test problem. The solution to this problem 
has a monotone decreasing density profile and therefore it does not display certain 
difficulties that may arise when the intermediate state has to be “built-up.” In 
figures 12-16 we present calculations for the Riemann problem 



(7.10b) (pi,q L ,P L ) = (0.445, 0.698, 3.528); (p Rl q R< P R ) = (0.5, 0, 0.571) 



used by Lax in [18]; see also [7], [9]. AH these calculations were performed with 100 
cells, h = 0.1, CFL = 0.8 and 85 time-steps using a component- wise reconstruction 
(5.11). In figure 12 we show the results of the “4-th order” ENO scheme using 
a component-wise reconstruction (5.11). Comparing these results to figure 11 we 
see that the component-wise reconstruction here is much “noisier” than in Sod’s 
problem. In figures 13 and 15 we show the results of the ENO schemes using 
characteristic reconstruction (6.3) for r = 2 and r = 4, respectively; comparing 
figure 15 to figure 12 we see that most of the “noise” in figure 12 is eliminated. 

In figures 14 and 16 we show the characteristic reconstruction R(x ; v n ) of the 
numerical solution in figures 13 and 15, respectively; R(x;v n ) is piecewise- linear in 
figure 14 (r = 2) and piecewise-cubic in figure 16 (r = 4). Tbe squares in these 
figures mark the values of R(x J+ , ± 0;v"); thus the difference between the two 
squares at the same location shows the size of the discontinuity in the reconstruction 
there (we recall that the circles in figures 13 and 15 are the values of fl(x_,;v n )). 
We see that the discontinuities in the reconstruction of the rarefaction wave are 
small enough to be graphically imperceptible. Surprisingly the discontinuities in the 
reconstruction of the contact-discontinuity are also rather small. Comparing figure 
16 to figure 14 we notice that the size of the discontinuities in the reconstruction 
for r = 4 is always considerably smaller than that for r = 2. It is interesting to 
note that even in the shock region in figure 16 (r = 4), the sum of the jumps in the 
reconstruction is only about 35% of the size of the shock, while about 65% of the 
shock jump is described by the smooth polynomial part of the reconstruction. 

We remark that because of the self-similar nature of the solution to the Rie- 
mann problem, the rate of convergence of any scheme is inherently limited to first 
order (see [27]). Comparing r = 4 with r = 2 in the solution of the above Riemann 
problems we notice a slight improvement in the smearing of the contact discontinuity 
(we have not used artificial compression in these calculations) and the description 
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of the rarefaction wave. Because of the self-similar nature of the solution it is better 
to compare the performances of two schemes by using x /t as the spacial variable 
and to find how many time-steps it takes to get well resolved intermediate states. 
Doing so for the problem (7.10b) we find that r = 4 with N = 35 gives about the 
same result as r = 2 with N = 70. 

(B2). Interaction of blast waves: In this subsection we present numerical ex- 
periments with the ENO schemes for the problem of two interacting blast waves: 



( u L 

0 

< x < 0.1 

(7.11a) 

u(l,0) = < IiM 

0.1 

< x < 0.9 


( ur 

0.9 

< x < 1 


where 

(7.11b) 

n = pm = pr = l. 9t . = l« = i« = o. P/. = io 3 , />« = lo- 2 , />« = io 2 ; 

the boundaries at i = 0 and x — 1 are solid walls. This problem was suggested 
by Woodward and Colella as a test problem; we refer the reader to [31] where a 
comprehensive comparison of the performance of various schemes for this problem 
is presented. 

In our calculations we divided the interval (0,1) into J cells by 
(7.12a) > = 1 '•' J ' 


where x ; marks the center of the ;-th cell. The boundary conditions of a so.id 
wall in x - 0 and x — 1 were treated by reflection, i.e. we defined auxiliary states 
v" , .... t/" r+ i for the left wall and v" +1 , ..., v n JJf . T for the right, wall by 
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We observe that representing the solid wall condition by the above reflection 
is very suitable for the characteristic reconstruction: V 3-wave approaching the 
right boundary is reflected as a 1-wave; consequently there is hardly any interaction 
\ between the waves in the characteristic variables (6.3a) and a situation of not having 

l 

| enough points of smoothness to choose from is thus avoided. 

\ In figures 17a- 17h we show the solution of the “4-th order” ENO scheme at 

f 

; t — 0.010, 0.016, 0.026, 0.028, 0.030, 0.032, 0.034, 0.038, respectively. We refer 

f rea d e r to figure 2 in [31] where a highly accurate solution is displayed and 

j a detailed description of the various interactions that occur at these instances is 

j' presented. The continuous line in figures 17a- 17h, 18 and 19 is the solution of 

the “4-th order” ENO scheme with J = 800 in (7.12a). Comparing this solution 
to the “exact” solution of Woodward and Colella in [31], we find that it shows 
all the important features of the various interactions and thus can be considered a 
“converged” solution. (The continuous line representing the solution with J = 800 is 
the piecewise-linear interpolation of {J*(xy; v n )}; consequently cusps in the solution, 
which do appear in R(x ; v"), are chopped in the graphic representation). The circles 
in figures 17a-17h show the values of R{ij\ v n ) of the “4-th order” ENO scheme with 
J = 400. Comparing the numerical solution for J = 400 to that of J = 800 we 
see that the velocity and pressure have already converged, while the density in 
figures 17g and 1 7h still deviates from the “converged” solution. This is due to the 
smearing of 3 contact discontinuities which are present in the solution at this time; 
the numerical results of Woodward and Colella demonstrate that the addition of 
contact-discontinuity steepeners” improve the density profile considerably. 

In figure 18 we show the solution of the “4-th order” ENO scheme with J = 200 
at the final time t = 0.038; In figure 19 we repeat the calculation in figure 18 for 
the “second-order” ENO scheme. Comparing figure 18 and figure 19 we see that 
the “4-th order” scheme gives a much better resolution. We remark that the results 
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Of the “4-th order” scheme with J = 100 (no, shown here) are of the same quality 
as those of the “second-order” scheme with J = 200. 

Wp nnf p n i sa 

■* w tiiicib d Ui 


paraooia interpolating 


/»(*_,) = 240, P(i„) = 0.01, />(*,) = 40 


has an mteiwal ,n whtch it is negative: the same is true for higher order interpolating 
polynomials tha, pass through these points. A situation of this type occurs in the calcu- 
late of the two interacting Hast waves just before the interaction in figure 17d, when 
the low pressure region in figure ,7c is shrinking to M computational cells. Since high 
order mterpolat.ng polynomial, may produce negative values of pressure and density in 
such drastic situations, we have imposed a -positivity condition” on the reconstruction 

step of our programs for the Euler equations. To ensure that R(x;v n , in the ;„h cell 
yields density and pressure that are positive, i.e. 

(7.13a) 

P > + E p y + V ti.\ (« - «,■)* 

k • dx k \ x=Xj k\ > 0 for |* - Xj\ < h/2 


fc=i 

we check whether 

r— 1 


(7.13b) 




5r r 


fc = i 

If condition (7., 3b) is no, satisfied we reduce the order of the reconstruction iocally 
at I z, until positivity is ensured. We observe that the LHS of the inequalities in 
(7.13b) is m in smooth regions, hence this positivity condition does no, reduce 
tbe asymptotic order of accuracy. Our computer program monitors occurances of 
order reduction due to the positivity condition; we have found tha, the order in 
the calculations of the -4-th order” ENO scheme has been reduced during two time 
Steps before the interaction in figure ,6d. and only at the interaction zone itself we 
have no, encountered any order reduction in the solution to the Riemann problem. 
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c Variants and extensions. 

c, CkaracUnsHc for <b -o/.r cose: In [15] « described an approx- 

imation to *(*,<), the solution-in-the-small of (4.2), which is obtained by trac.ng 


ipproximate chaiacteristics to the 


initial data. This approximation 5(*,t) can be 


extended to arbitrary order of accuracy 


follows: Let a" A denote 


(7.14a) a",. -°iO. R(* 3+i +0-,v")) 

where S (u, , « 2 ) is defined in (4.7), and let a(x) denote the interpolation of ^ j by 
H m (3.1) with m = r - 1, i.e. 


(7.14b) 




(7Hc) 5(x) = m = r-l. 

The approximation e(x,t) is obtained by prescribing constancy of the solution 
along the approximate characteristic lines 


(7.15a) 


x = Xq + fl(*o)^ 


(7.15b) 


v(i 0 + a(x o )t,0 = ^C x o,o) = fl(zo;0; 


(7.15c) v{x t t) = RMx,t)X) 

where *„(*,<) is the solution of the algebraic equation (7.15a). Let denote 

the solution to (7.15a) for m = l in (7.14c); if « and < are such that 


(7.16a) 


ly_x(0 < Z < *>+ x(0 
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where 


(7.16b) 


then 


(7.16c) 


*»+*(0 — I i + a + t 


*o(*.0 = j 


x i±j ~ x J~k 


(X - *y_l(<)). 


x i+ x(0 -Xj_i(0 

For m > 1 we obtain x 0 (x,f) by solving (7.15a) with Newton- Raphson iterations 
starting with the initial guess (7.16c). 

Using v(x,f) (7.15c) we define the following variant of (4.14): 


(7.17a) 


V? + 1 = Vy — A(/y + l - /j_l) 


(7.17b) 


/y+* = £ afc /( 5 (*;+i>& T ))- 


*=0 


We have started the development of the ENO schemes with the version (7.17); 
later on we have replaced the characteristic method by the Cauchy-Kowalewski 
procedure which offers a unified approach in extending the scheme to include forcing 
terms and to systems of conservation laws. Our numerical experiments show that 
the two versions are computationally equivalent, although the version with the 
characteristic method (7.17) seems to be slightly more accurate than (4.14). 

We remark that the scheme (7.17), as the scheme (4.14) with f ROE in (4.21a), 
also admits any discontinuity with a(u£,,u /?) = 0 as a stationary solution. This 
can be easily rectified by replacing the “shock curve” x l+ i(0 in (7.16b) by an 
appropriate fan. 

C2. Semi-discrete formulation and Runge-Kutta methods: The semi-discrete 
version of the ENO schemes can be derived either directly from (1.4) or by letting 
r — * 0 in (4.14), (5.9). It takes the form 


(7.18a) 


f t M0= 4[/,>j(0-/>-*(0|s<j *,(0 
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where 

(7.18b) 4 +x(0 = / i? (/2(i y+ i -0;v(<)), £(*y + x + 0;v(f))); 

here Vj(t ) is an approximation to u(xy,f);v(f) = {vy(f)}; /^(«i,u 2 ) is either the 
exact flux (4.5) or f ROE (4.21). (5.11). 

Considering (7.18) to be a system of ordinary differential equations in t for the 
vector v(t) = {vy(f)}, we can solve the problem by using a numerical ODE solver. 
In [2] we present two sets of numerical experiments in which we use Runge-Kutta 
methods of appropriate order to approximate the solution of (7.18). In the first 
set of experiments we apply the scheme to the Riemann problem (7.10a) for r = 1, 
2, 3, 4, 5, 6. In the second set of experiments we apply the scheme with r = 2,4 
to a Laval nozzle problem which involves the addition of a forcing term to the 
Euler equations (6.1). In these calculations we have used RP, f ROE and CFL = 
0.5. Comparing the results of the Riemann problem to those in the present paper 
we find them to be of similar quality. The numerical experiments of (2| indicate 
that the semi-discrete formulation (7.18) with Runge-Kutta temporal discretization 
does not generate spurious oscillations for CFL < 0.5; however when we increase 
the CFL number beyond 0.5 we start getting some oscillations and eventually the 
scheme becomes unstable. 

The main advantage of using the Runge-Kutta temporal discretization is the 
ease of its programming; however it seems to be less efficient than the fully discrete 
formulation and also requires more storage. 

C3. Variable grid and front tracking: In section 3 we have pointed out that the 
non-oscillatory interpolation H m (3.1) - (3.5) and the reconstruction via primitive 
function (RP) (3.6) - (3.10) are well defined for non-uniform grids, see appendix. 
Since the solution-in-the-small step also does not require uniformity of the grid, we 
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may 

by 


compute new cell-averages v? +I in (1.14c) on any choice of intervals (J? + '} 


(7.19) 


v"*' = — L- f fk (i,(„+ 1 -0)<fi; 
1 |l“ +1 |7r-+‘ 


here /' = (f‘_ r, «' + a) and |Jj| = €' +i -<£-*• Wm « lhe same rali ° nlU “ bef ° rC ' 
our approximation to (7.19) becomes 


(7.20a) 


|J7 +, |v“ + ‘ = |f7|v”-r(/J +i -^ 5 )- 


The numerical flux /f i is consistent with /( u) - Uj+s“ 

J“T 3 


(7.20b) 


and can be expressed as 






(7.20c) = 52 <*kf R {V]{ X j+$ d /5fc^j+A.^fc T )i W J + l( I J+i +/ 5 fc T<T i4 .L,PkT),<?j+L) 

5 fc =0 


where 


(7.20d) / R (tti,tta;<r) =/(V r (a; u i, u a)) ^^(a; Ui, u 2 ), 

we recall that V(<r; up, »„) denotes the value of the solution to the Riemann problem 
(4.3) at x/t = a. Roe’s linearization (5.11) yields the following approximation 

(7.21) / r<,£ (u,,u 21 <t)= j|/(u,) + /(»j) -u(“i +“s) 

rn 

- 52 5fc( u i, u 2)i< i fc(^) - <T i r fc(*)i- 

<c=i 

In figure 20 we show the results of the scheme (7.20) with j R approximated 
ky jroe (7.21) for the Riemann problem (7.10b); the values of (£y } in this 

calculation were chosen by the self-adjusting grid algorithm of (13). This algorithm 
provides an automatic way to place interval end-points £" + l at the location of 
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significant discontinuities and thus avoid their smearing by the cell-averaging step 
(7.19). The calculation in figure 20 was initialized by taking the exact solution of 
the Riemann problem at t = 0.5 (at which time there are 4 grid points between 
the contact-discontinuity and the shock). The results displayed in figure 20 show 
the numerical solution of the scheme with r = 4 after 100 time steps with CFL 
= 0.5. These results clearly demonstrate the adaptability of the ENO schemes to 

front tracking techniques. 

We note that the use of irregular grids disallows the extra order of accuracy 
which was gained in ( 1. 18) for a uniform grid. Numerical experiments with irregular 
grids (where «” +l is randomly selected within a specified interval) show that the 
error in solving the scalar smooth problem (7.1) by the scheme (7.20) is O(V') in 
the I, , L , and L„ norms. However comparing figure 20 with fig-are 15 we observe a 
considerable gain in resolution in spite of the reduction in formal order of accuracy. 

C4. Extension to 2D: In this subsection we outline the extension of the ENO 
schemes to the solution of the two-dimensional IVP 

(7.22) u* + /(«). + ?(«)» = °> u ( z > °) = Uo(l ’ y) - 


We note that Strang-type dimensional splitting [29] is only second-order accurate in 
time, and therefore is unsuitable for extending the higher order accurate members 
of the ENO schemes to 2D. 

Let ti ) denote the two-dimensional “sliding average of tv 


(7.23) 


. <• Ax/2 /• Ay /2 

^( I)V ) = — — - / w{x + + l) dr l d Z’ 

{ AxAy J-ax/2 J- Ay /2 


Integrating (7.22) over the computational cell I tJ x [< n ,<„ + ij, /,j - 




j+ji’ we 


find that = S(x„ *n) satisfies the equation 


(7.24a) 


- A .(/,+ *.,- - /.-*.,) - A ,(fc.y + i - 9i.j- 


■*) 




I 
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where A x = r/Ar, A y = r/A y and 

- 1 f T f&v / 2 

(7.24b) /»+*, i = ^ J-&y/2 ^ Xi+ 2 ,y * + + *))*?<**» 

j /»r /•Ax/2 

(7 - 24c > = 7 a7 / / »(«(* + <. %+*.». + <))<*«*• 

mz Jo J- Ax /2 3 

The abstract form of the ENO schemes for the solution of (7.22) remains (1.10), 
i.e. 

( 7 -25) v" +1 = A k E(t) ■ R(; •; v n ), v° = S 0 . 

As before E(t) is the exact evolution operator of (7.22); however, A h is now the 
2-dimensional cell-averaging (7.23) and R(x, y; tf) is an appropriate ^-dimension J 
reconstruction of tu(z,y). In the scalar constant coefficient case 

(7.26a) u t + au x 4- bu y - 0, u(z, y, 0) = u 0 (x, y), 

the ENO scheme (7.25) becomes 

(7.26b) v" +1 = R(x x - ar, y y - 6r; v n ), v? = So(z„ y y ). 

In [12] we present numerical experiments with the ENO scheme (7.26) for 
the scalar constant coefficient case, where the reconstruction R(x , y; tZ>) is obtained 
via a two-dimensional deconvolution. Expanding w(x 4- £,y + rj) in (7.23) around 
£ = rj = 0 we get as in (3.12) 

(7.27) t5(z, y) = w{x, y) 4- a 2 [(Az) 2 u/ XI + (Ay) 2 w ys ,j 4- ^[(Ai) 4 ^^ 

+ 2(Az) 2 (Ay) 2 u; xiyy -f (Ay) 4 tu yyyy j 4- 0(A 6 ) 

Multiplying both sides of (7.27) by (Az) fc (Ay)'~ fc and truncating the ex- 

pansion in the RHS at 0(A r ), we get as in (3.13) an invertable system of linear 
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equations which expresses « and its derviatives in terms of w and its derivatives. 
We set 


(7.28a) 


{D 00 )ij = w{xi,yj ) 


and obtain approximations 


(7.28b) (£*•%= (A*)*(Aj,)' ak *‘ . 


fckgj «(*(,!/,) + 0(A P ), 1 <* + t<r-l; 


then, as in (3.17), we invert the system of linear equations to get the following 
approximations to w and its derivatives 


(7.29a) (Z>*% = (Ai )*(A s,)'-L_ W ) + 0(A P ), 


dx k dy l 

Using (7.29a) we define A in the cell by 

r-l . / 


0<fc + /<r — 1. 


(7.29b) S(x, y ;t5) = VlV ( D‘'-‘ )l ,Y'U£z£l)Vs'-y J V-‘ , . . 

a * ) \~Kf) • (*> »)£/*. 


The approximations D* ‘ in (7.28b) are obtained by a sequence of applications 
of the one-dimensional operation (3.15b), which we rewrite now in the following 


operator form: 


(7.30) 


{G ' m * U)j ~ ~ ° : u > - 7n H ^j + 0; u)); 


dz l 


her, « denotes the one-dimensional vector {*(*,)}', , . Using (7.30) and the notation 
convention = (tS(i,, v>)} , we define 


(7.31a) 


{D k '°) i j _ (A x) fc (G* • t5.. y ),-, 1 < Jb < r - 1, 


(7.31b) 


( D °% - (&y)‘( G 'm • «"..•) y , 1 < l < r- 1 . 
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To obtain approximations to the mixed derivatives of w we first evaluate 

(Df-'ji, = (Ari'lGj,,-* • 1 < I < r - 1 - k, 

(£>!;■% = (**)*[<&_, • 1 < fc < r - 1 - i, 

and then define 

(7.31c) {D k l )ij =M{{D k l l ) i j, ( D kl )ij ), 


where M is the min mod function (3.16). 

We observe that the restriction of the two-dimensional reconstruction (7.29) 
to y = j/j, i.e. R(x,yj ;& ) is identical to the one-dimensional reconstruction (3.18) 
applied to the restriction of tiD to y = y } f i.e. R(x\ <D(«, yy)); the same observation 
applies to the restrictions to x = x,-. 

We recall that the one- dimensional reconstruction is essentially non-oscillatory 
only if discontinuities are seperated by at least r + 1 points of smoothness. In 
the one-dimensional system case we had to overcome the problem of collision (in 
time) of discontinuities; in the two-dimensional case we also have to worry about 
intersections (in space) of curves of discontinuity. In order to study the severity 
of the problem we have experimented with the constant coefficient problem (7.26a) 
with the initial data 


(7.32) 



(*.y) € 5; 

( x,y)eU-S 


here U = [—1,1] x [—1,1] and 5 is a rotated square contained in U . In [12] we 
present numerical results which are obtained by applying the scheme (7.26) with 
r = 1,2,3, 4 to the initial data (7.32) with periodic boundary conditions on dU . 
These results show that indeed small spurious oscillations are generated for r > 2 
at the corners of 5’, however it seems to us that they are small enough to be 
computationally acceptable. 


M 
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then we evaluate 


(A.7) 


It is easy to see that 


do Z = 1> r 

Su = Si-u-iZi-i 

do k = Z + 1, t 

St.k = Si.k-I + Si-i,k-iZk-i 

r 

q^(x e ) = 


fA8 i = 

V ' ’ l=k 

We note that the algorithm (A.6) • (A.8) is defined for a non-uniform grid. 

When the grid is unifora we can obtain (A.8) in two steps: First we take *„ = Vi 
and observe that Z, = -Ih in (A.6); consequently {$,*} are independent of .. 

Denoting 

(A . 9) i k = h k ii.„, c fc = 

and using the convention 1„ = 0 for k > r, we get for 1 < r < 6 

' c 0 = u(y<) _ 

Ci = d, - 4/ 2 4- 4/ 3 - 4/ 4 4- 4/5 -4/6 

C3 = d 2 -d z - i- 114/12 - 54/6 4- 1374/180 
( A 10 ) < c 3 = 4 ~ 1-54 + 1-754 - 1-8754 

c 4 = 4 — 24 + 174/6 
c 5 = 4 - 2.54 
< Co = 4 


«<*) = E i 


c k /x- y» 


— t/_- \ k 


Thus 


(A. 11) 


Jt! V h 


WM - 1 (TTji 


I — r* 

fiecc.nsiruct.on « /unci, on ( RP In this case /, = is the 

j-th cell and I, = t, = *(»/ + »s>.) » « ntCT - The « iveD daU * 


( A. 1 2) 
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from which we evaluate the point values of the primitive function 

fc-i 

(A. 13) w{y k ) = ^(l/ 1 +i - yi)m . 

1 = 0 

Applying the algorithm (A.6) - (A.8) to 

(A. 14) di, k ~W[yi t ...,yi+k\ 

with i = i(j) selected by (3.4), we obtain the values of <7 (/) (z y). Using the definition 
(3.8) in (A. 4) we get the coefficients of the Taylor expansion in (A.l) by 

(A. 15) = 

We note that when the grid is uniform y k = ifc- 1/2 we caQ a l so use 
algorithm (A.9) - (A.ll). 

Reconstruction via deconvolution (RD): We recall that RD is used with a uni- 
form grid so that the given data t Dj can be thought of as point values tn(zj) of the 
sliding average function (1.3). Applying the algorithm (A.9) - (A.ll) to 

(A. 16) d,. fc = u >(z.', .... z«+fc] 

with x c = Xj we get in (A.ll) for » = i(j - 1) the values of 

d k 

(A. 17a) h k -^ H r (xj - 0]w)\ 

when we apply this algorithm with i — i(j) and x,. = xy we get in (A.ll) the values 
of 

d k 

(A. 17b) h k H r {xj + 0;u>). 

Next we evaluate Z> fc .y in (3.15) by taking the min mod of the appropriate values 
in (A. 17a) and (A. 17b). Finally we use the back-substitution (3.19) to obtain the 
coefficients of the Taylor expansion (A.l) 

h* = Jj D k .,lh k . 


• •*- -ix .-JT.-Jgrc 



(A. 18) 
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We remark that the use of the algorithms (A.9) - (A.ll) is preferable to that 
of (A.6) - (A. 8) since it enables us to save computing time by rearranging the 
operations (A.16) - (A.17) as follows: First we set i = i{j) in (A.16) and evaluate 
(A.9) - (A. 10). Using the same coefficients c k in (A.10) we now apply (A.ll) to 
x c = Xj and x c = * i+1 to obtain H r (x j + 0; u>) and £* H r (x j+l - 0;®), 
respectively; the min mod operation (3.15) is then performed in a following sweep. 



1 


J 

s 
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Solution of the periodic IVP (7.1) at t = 0-3 by ENO 
schemes based on RP . 
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Table 2a. Solution of 












Table 2b. Solution of 
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Figure la: Second-order solution (t - 0.3) 



Figure 2a: Second-order solution (t » 2/ff) 


Figure lb: Fourth-order solution (t = 0.3) 



Figure 2b: Fourth-order solution (t « 2/*) 










Figure 6®: First-order solution with 
' E (. J = 40, N = 80) 



Figure 6b: Second-order solution 

th f ROE ( J = 40, N = 80) 


Figure 6d: First-order solution with 
fROE ( j - so, N = 160) 



Figure 6e: Second-order solution with 

fROE ( j = go, S = 160) 



i ** * , 

Figure 6c: Fourth-order solution with 

RaE ( J = 40 , iV = 80 ) 





Figure 7«: /( u) (7.7b) and its con- 
vex hull for u L = -3, u a = 3 


Figure 7c: Second-order solution with 
exact f R ( J = 40, N - 20) 
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Figure Tb: First-order solution with 
exact f R ( J = 40, N = 20) 



Figure 7d: Fourth-order solution with 
exact f R (J = 40. <V = 20) 
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Figure 8a: First-order solution with Figure 8b: Second-order solution 

B ( J = 40, N = 20) with f ROE {J = 40, N = 20) 



Figure 8c: Fourth-order solution with 
(J = 40, = 20) 
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Figure 9s "Second-order" solution of Sod’s prob- 
lem (characteristic reconstruction). 
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Figure 10: “Fourth-order" solution of Sod’s prob- 

lem (characteristic reconstruction). 
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Figure 11 
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“Fourth-order" solution with component- Figure 12: “Fi irth-order* solution with component- 

w.se reconstruction of Sod’s problem. «i» reconstruction of Lax’s problem. 
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■Second-order* solution of Lax’s prob- 
lem (characteristic reconstruction). 



Figure 14: Characteristic reconstruction of the 

■second-order* solution in Figure 13. 
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Figur. 17®: 1 = 0.010 Pig»r. 17b: i = 0.016 
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Figure 17: “Fourth-order 7 ’ solution of the interacting blast waves problem; J 

continuous line; J = 400 - circles. 
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Figure 20: “Fourth-order* solution of Lax’s problem with a self-adjusting grid 
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